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SECTION I: INTRODUCTION 



Historical Notes 



To most engineers practicing in industry the high speed digital computer 
is an unknown quantity. Just five years ago an engineer's chances of having 
an introduction to computing while in college were slim indeed. 

Today virtually all the major engineering schools in this country have access 
to digital computers. Within the next five years the engineer who graduates 
without an introductive course in the use of computers will be the exception. 

The growth in the use of the digital computer has been tremendous . The 
first large scale digital computers appeared in the early 1950s as an out- 
growth of interest generated through the use of punched card calculators 
such as the IBM Card Programmed Calculator. Over 4,000 digital com- 
puters were installed in I960, varying from small engineering machines 
to very large commercial and scientific data processing systems. 

Like all technological developments, the digital computer has its historical 
antecedents. Computing itself is one of the oldest human activities, re- 
quired by all civilizations for the conduct of business and the development 
of sciences. One of the most practical inventors of recent times, Charles 
Kettering, inventor of the auto self-starter, made the statement, "You 
cannot understand a thing unless you can give it a number." In many ways 
civilization has always been tied to this need for quantitative or numerical 
description. 

In fact, the oldest surviving written document is a set of business records 
kept by a Sumerian of Mesopotamia 5, 000 years ago. 

The development of computing methods was painfully slow. Perhaps this 
was due in part to the unfortunate choice of number symbols (e.g. , Roman 
numerals) and of number bases (e.g. , the Babylonian choice of the base 
60). Gradually the Arabic numeral system and the base 10 won out as 
modes of computation, although vestiges of other systems remain impor- 
tant (e.g. , Roman numerals for date inscription, and base 12 for the inches 
in a foot measure). 

Perhaps, too, this slow march forward in computing methods was due to 
the fact that they were not greatly needed, for certainly the first big 
breakthrough occurred when the science of astronomy depended on tremen- 
dous computing efforts for its advance. In the 16th century these efforts 
were greatly eased by the development of logarithms by Napier and Briggs. 
Considering the size of numbers used in astronomy, one can easily see why 
the logarithmic transformation which allows multiplication and division to 
be handled as addition or subtraction was a significant contribution. 

In 1642 Blaise Pascal, a French mathematician, built the first successful 
digital computing machine . While this device demonstrated the feasibility 
of mechanical computation, the use of such devices remained very limited. 



A most interesting personality in the science of machine computation was 
the misunderstood genius Charles Babbage of England, whose major work 
was done in the middle 1800s. During his lifetime, Babbage conceived 
and designed on paper what he termed the analytical machine. In virtually 
all respects this machine was to perform like the modern digital computer, 
except that it was to be mechanical rather than electronic . The fasci- 
nation which Babbage's work holds for present-day computer people can 
be explained when some of his ideas are discussed. 

For example, his analytical machine was to have the ability to "remember" 
1, 000 numbers up to 50 digits in length. This capacity is greater than 
that of many electronic computers now in existence. His machine was to 
be controlled by punched cards, an idea which he borrowed from the use 
of cards in controlling looms. In the light of the present day role of the 
punched card in data processing, Charles Babbage seems a prophet. 

It is a certainty that if Babbage had been content to build on a smaller 
scale rather than design on such a grandiose level, the science of com- 
puting would have profited greatly in the 100 years between his analytical 
machine and the present day. As it was, however, Babbage never com- 
pleted his machine, and in general we are aware of his genius primarily 
through the writings of an admirer, one Lady Lovelace, who wrote most 
descriptively of the design and philosophy of the analytical machine. 

In truth one of her statements might stand as a preface to the study of the 
modern electronic computer: "The Analytical Engine has no pretensions 
whatever to originate anything; it can do whatever we know how to order 
it to perform. " 

Charles Babbage did successfully build what is termed a difference machine 
for the computation of tables. Many of the early efforts in machine com- 
putation took this line of construction of machines to do particular comput- 
ing jobs. A simple example will illustrate a principle made use of in such 
machines: 

To construct a table of values for 

F = M 3 + M 2 

where M is to take on the values 0, 1, 2, 3, , set down the first 

few entries and take successive differences in F: 



Differences 



M F 1st 2nd 3rd 









1 


2 


2 


12 


3 


36 


4 


80 


5 


150 


6 


252 



2 


8 




10 


14 


6 


24 


20 


6 


44 


26 


6 


70 


32 


6 


102 







The third difference is a constant. The table can be completed by working 
back from the third difference to compute further entries of F by simple 
additions . It is relatively simple to build a machine to perform these 
simple additions to successive columns. 

The principle used here has been to take successive derivatives of the 
function F. 

dF 9 

d 2 F 

z-^ = 6M + 2 

dM 2 

d 3 F a 

o = 6 

dM 3 

The third derivative is the constant 6. 

A great many of the functions for which tables are desired in the fields of 
navigation, astronomy and surveying can be handled by such techniques. 
The literature of computing abounds with the application of difference 
procedures. 

However, modern computers are a far cry from the difference machines 
of Babbage's day. 

At about the same time Babbage was laboring with his "Analytical Engine, " 
George Boole, an English mathematician, was laying the foundation for 
later developments with his study of the algebra of logic. Boolean algebra 
is the foundation of work in machine logic (so necessary in the design of 
digital computers) and has contributed in the logic of business contracts 
and law (equally necessary in many areas of data processing). 

In the 1930s the work of modern mathematicians, in particular the late 
John von Neumann, spurred interest in the possibilities of machine com- 
putation. Various computers — both special-purpose analog devices and 
true general-purpose electromechanical digital computers — were built 



at universities. Harvard and MIT were the centers of most of this activity. 
World War II was a catalyst to machine computing in that the military 
technological advancements demanded tremendous computational facilities. 
By the end of the war the logic for construction of a general -purpose 
digital computer was fairly well specified. 

As makers of computing devices for accounting and business, manufac- 
turers of business machines became involved in this effort of machine 
computation. The first large scale general -purpose digital computer was 
built by IBM for use in scientific research. Known as the IBM Automatic 
Sequence Controlled Calculator , it was presented to Harvard in August 
1944. Its capacities were not surpassed until the IBM Selective Sequence 
Electronic Calculator was dedicated in January 1948. 

The SSEC was dismantled to make room for its successor, the IBM 701, 
announced in March 1953. At this time the general belief was that a few 
large scale machines could handle all the computing needs of the country. 
The demand for more and more machines when the first few became 
available soon made apparent the vast needs for computation in both tech- 
nology and business. Today these needs have been multiplied by further 
recognition of existing needs and the discovery of new approaches to old 
problems. 

How does this affect the engineer? 

The field of engineering analysis for computers has experienced tremen- 
dous growth in the past five years. The advent of the electronic computer 
has created a new approach to engineering problems. The costly "build 
and try" method is now often replaced by "construct mathematical model 
and simulate." 

In some respects this represents a return to the textbook for the engineer. 
He must describe the physical problem by means of a mathematical model. 
Then the computer can operate the model and describe operating results. 
The real power of the computer is that it will quickly describe operating 
results for many sets of operating conditions. 

Many articles have been printed in professional publications on the solu- 
tion of difficult engineering problems by the use of electronic computers . 
These applications have only scratched the surface of possible engineering 
use, and a lack of knowledge and experience prevents the full use of the 
computer as an engineering tool. The engineer is simply not aware of 
what a computer can do for him . It is highly desirable , therefore , that 
every engineer understand how problems can be described for handling 
by electronic computers. 

It is convenient to think of computer problems as falling into one of two 
classes: 

1. Obvious computations 

2. Iterative problems 



Obvious computations are those which have in the past been handled by- 
slide rule or paper and pencil. They also include those for which the tech- 
nique of solution, while known, has required too much computation to tackle. 

Iterative problems comprise a large area, mostly unexplored. To explain 
this, the question should be asked, "What is the function of the engineer 
today?" 

The engineer is a person who builds physical systems to do particular 
jobs. If a particular construction does not work, modifications are made 
until it does work. These modifications are often the result of vague 
ideas. Engineering is largely by trial and error. 

A very obvious conclusion which one reaches in talking with engineers is 
that they employ a minimum of their textbook analysis techniques . Math- 
ematical approaches soon pick up dust on the engineer's worktable. Why? 
One reason is that real life simply does not fit the pat equations of the 
textbooks. Often the basic difference is that in a textbook it always seems 
possible to apply mathematical techniques and arrive at solutions of the 
kind 

x = f (y, z) 

where the value of the variable of interest can be determined simply by 
plugging in known values for y and z on the right-hand side. However, in 
real-life engineering, often the problem becomes stated as 

x = f (x, y, z) 

where x is on both sides of the equation, x cannot be solved for explicitly; 
that is, there is not enough information to remove x from the right-hand 
side. The problem is too complex. 

This problem can usually be handled mathematically by estimating the 
value of x, testing the estimate, and, if it is wrong, making a better 
estimate. Essentially, numerical analysis (in computer mathematics) is 
the science of making progressively better estimates. To estimate 
repeatedly is to "iterate." This iterative idea is the heart of the computer 
approach, whether for the simplest cam problem requiring algebra, or 
for the most complex nuclear problem involving sets of differential equa- 
tions . 

This manual describes typical problems which have been handled by engi- 
neers making use of computers . No attempt is made here to give all the 
mathematical techniques required in computer analysis , nor to explain 
fully the programming techniques required for control of electronic com- 
puters. It is felt that an introduction to actual computer problems, with 
an attempt toward classification as to type, will give the practicing 
engineers insight into what computing can do and perhaps suggest possible 
uses they can personally make of engineering computing analysis . 



The Digital Computer 



The layman's image of the electronic computer contains a pinch of mystery, 
a dash of magic and a number of false analogies to human beings. These 
misconceptions have in some cases been nurtured by computer manufac- 
turers through their use of such terminology as 

1. Memory — for that part of the computer which stores information. 

2. Nerve Center — for that part of the computer which controls its 
operation. 

3. Brain — for that part of the computer which accomplishes arithmetic 
and logical operations. 

While these analogies do help in understanding the functions of some of the 
components of a computer, they also tend to develop a camouflage of 
complexity . 

The computer, like the automobile, airplane, nuclear reactor or any other 
engineering achievement, is the product of a group of engineers. (It should 
not be surprising that in developing this product, engineers have used other 
computers to solve some of their design problems . ) As in the case of the 
automobile and airplane, a person need not be capable of designing and 
building a computer in order to use one. Few persons who drive an auto- 
mobile are capable of describing the complex of gears , rods and shafts 
that allow a rotation of the steering wheel to change the direction of travel. 
This, however, does not prevent them from becoming excellent drivers. 
Similarly, a general knowledge of the functional components of a computer 
is all that is needed to use one effectively. Once the basic concepts are 
understood, greater familiarity is easily developed through experience. 

What are these basic concepts? To answer this, consider the basic com- 
ponents of a computer. 

STORAGE 

That part of a computer which most differentiates it from a calculator 
(desk-type, slide rule, adding machine) is its storage — or, to use the 
unapproved analogy, its "memory." 



STORAGE 



A computer can store numbers, alphabetic characters and some special 
symbols. In any calculation, it is necessary to store the numbers which 
begin the calculations, intermediate results and — at least temporarily — 
final results. A desk calculator has a very limited ability to store data. 
The operator usually must enter each number as it is needed. On the 
other hand, a small computer (such as the IBM 1620) can store 20,000 



digits and call them out of storage whenever they are needed in a calcula- 
tion. 

THE STORED PROGRAM 

The ability to store large amounts of data that is to be used in or has been 
developed by calculations is only part of the function of storage. Equally 
significant is the fact that the computer also contains within storage the 
program for the calculations to be performed. 



Data 

STORAGE 

Program 



A program is made up of instructions which, when acted upon by the com- 
puter, cause it to go through the sequence of operations necessary to arrive 
at a meaningful result. A single instruction may: 

1 . Cause data to be brought into storage from some external source , 
such as a card reader. 

2. Cause a specified arithmetic operation to be performed on selected 
numbers . 

3. Constitute a logical test to determine what part of the program should 
be performed next. 

4. Cause results to be sent from storage to a recording device, such as 
a typewriter. 

After both the data to be operated upon and the program which describes 
the operations are in storage, the computer is free to proceed with a 
series of calculations at its own natural speed. For present-day machines, 
this may be from 50 to 500,000 additions per second. 

This leads to another salient fact: A computer can change or modify its 
own program! Since the program is stored in much the same form as 
data, one portion of a program may be a sequence of operations which 
will examine another portion of the same program and change it by arith- 
metic or logical manipulations. What advantage is to be gained from this? 

First, it means that one set of instructions can be used to operate on a 
number of sets of data stored in different locations in storage . As each 
set of data is processed, the program is changed to refer to the next set 
of data. 

Second, it allows the program to be changed on the basis of conditions 
which arise during the calculations. For example, suppose a calculation 
involves the evaluation of 



y =x< 



, for x ^ 1 



y = -x +4x-2, for x > 1. 

At the point in the calculations where the value of x becomes greater than 
1, that portion of the program involved in the evaluation of y can be changed 
to use the second of the two equations . This change might be accomplished 
by replacing one set of instructions with a new set stored elsewhere for 
that purpose, or it could also be done by causing the computer to select 
one of two instruction sequences on the basis of whether x is greater than 
1 or not greater than 1. 

ADDRESSES 

To be able to select the item of data or the instruction to be operated upon 
next, the computer must be provided with a means of locating the desired 
information in storage. For this reason, storage is divided into units and 
each unit is identified by an address. Different computers are designed 
with different-size units of storage . The basic unit may vary from one 
decimal digit (as in the IBM 1620) to 36 binary digits (as in the IBM 7090). 
In other computers, units of one alphameric character (that is, a decimal 
digit, letter of the alphabet, or special character) or ten decimal digits 
are used. 

Whatever the size of the basic unit, each is assigned a numerical address 
which identifies the location. Manipulation of the data stored in a location 
is accomplished through the use of the address corresponding to the loca- 
tion. 

CONTROL UNIT 

The control unit of a computer provides the means whereby the stored 
program causes the desired operations to be performed in the sequence 
specified. The control unit reads an instruction from storage, examines 
it, sets up the circuit conditions to perform the operation, and, when the 
operation is completed, reads the next instruction from storage and re- 
peats these steps. 



CONTROL 



STORAGE 



Generally, the first operation to be performed is manually entered into 
the control unit by the machine operator. Thereafter the action of the 
computer is completely controlled by the program in storage as inter- 
preted by the control unit. An instruction which tells the computer to 
stop can be the last one in the program. 



ARITHMETIC UNIT 

This component contains the circuitry which performs arithmetic on num- 
bers taken from storage. It usually includes a limited amount of storage 
in which to hold the operands involved in the arithmetic . 



CONTROL 




At this point it is wise to pin down the kind of arithmetic accomplished in 
a computer. The two basic types of computers — analog and digital — 
have different ways of accomplishing arithmetic. The analog computer 
performs arithmetic by measuring, the digital computer by applying a set 
of rules to the numbers involved. 

To illustrate , suppose that two numbers to be added are 56 . 8 and 29.5. 
A simple analog device for obtaining the sum would be a ruler divided 
into 100 units. The ruler is laid on a piece of paper and marks are made 
on the paper next to the and 56.8 points on the rule. Then the ruler is 
picked up and the zero point of the ruler is placed next to the mark on the 
paper representing the 56.8 point. Now a third mark is made on the paper 
opposite the 29.5 point on the rule. All three marks, of course, are made 
to fall on a straight line. Finally, the distance between the first mark and 
the third is measured with the rule . Because of the possibility of measure- 
ment errors inherent in such a procedure, the result obtained might per- 
haps be 86.2. 



56.8 



J I I I I I I L 



n r 



T 1 1 1 r 



29.5 

X 

1 1 r 



86.2 



A schoolboy will add the numbers by writing one under the other and ap- 
plying a set of rules involving such things as the sum of two digits and how 
to handle the carry when the sum is greater than 9. Thus, 

11 
5 6.8 
2 9.5 
8 6.3 

This is the technique used by the digital computer — that is, the digital 
computer performs arithmetic by the application of a set of rules. There 
is no error inherent in such a technique. 

BINARY ARITHMETIC 

While on the subject of arithmetic, it is well to dispense with another 
question which is bound to arise in the reader's mind: What is binary 
arithmetic ? 

A problem in computer design that the engineer may have recognized by 
now is that to store and manipulate decimal digits requires some device 
capable of assuming ten stable and unique states — one state to represent 
each of the digits to 9. While such devices are available (the notched 
wheels in a desk calculator are an example), they lack the speed of opera- 
tion and small size necessary for a computer. Furthermore, they don't 
fit into an electronic system. 

Many electronic devices are available, however, which can assume two 
stable and unique states. A tube, for example, can be conducting or non- 
conducting; a magnetic field can have clockwise or counterclockwise 
rotation; a pulse can be transmitted or its transmission can be blocked. 

The computer designer has at his disposal, therefore, a number of devices 
and techniques for operating in a number system with the base 2 — the 
binary number system. The rules for arithmetic are also much simpler 
in the binary number system than in the decimal system (base 10) which 
we are accustomed to using. 

ADDITION TABLES MULTIPLICATION TABLES 
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10 
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1 



As a result, there are two basic types of digital computers: one which 
operates entirely in the binary number system and another which codes 
decimal digits as binary numbers. For the latter, the decimal digits 
appear as their binary equivalents: 
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DECIMAL BINARY 



DECIMAL BINARY 









5 


101 


1 


1 


6 


110 


2 


10 


7 


111 


3 


11 


8 


1000 


4 


100 


9 


1001 



As far as the user is concerned, it is not necessary to be familiar with 
the binary nature of the computer. The computer is capable of translating 
the engineer's language to its own before starting a program and performs 
the reverse operation when communicating its results to the engineer. 

INPUT AND OUTPUT 

So far, the methods for getting data and programs into the computer and 
for getting the results out have been ignored. Every computer must have 
a means for communicating with the outside world. 

Typical input devices are punched card readers, paper tape readers and 
manual keyboards. Card and paper tape punches, typewriters and printers 
are typical output devices . Magnetic tapes and magnetic disks provide a 
means of storing data 



CONTROL 
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STORAGE 






INPUT 


OUTPUT 
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rn 



ARITHMETIC 



externally from the computer in a form which allows rapid re-entry. 

The use of input-output devices is under control of the stored program. 
For example, an instruction to read a card will cause the card reader to 
start up, feed and read one card and transmit what has been read into 
storage . 

SUMMARY 

Thus there are five basic components of a computer: 

1. A place to store data and instructions in such a way that each item can 
be located when needed. 
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Exercises 



2 . A means of controlling the overall action of the system . 

3. Devices for performing arithmetic and other operations required 
in manipulating the data to obtain desired results. 

4 . A way of getting information into the computer . 

5. A way of getting results out of the computer. 



1. The decimal system has as its base the number 



2. The numeral system commonly used is the 

numeral system. 

3 . Logarithms aid greatly in doing the arithmetic operations of 

and when working with many significant 

figures . 

4 . Difference principles were used in the construction of 



5. The first large scale computers were electromechanical in design; 

the middle 1950s saw the common use of the design, 

and the computers being designed and built now have solid-state 
(e.g., transistorized) components. 

6 . The general iterative problems may be stated by a formula of the 
type 

7 . Three computing devices with which most engineers are familiar 
are , , and 



8 . One of the difficulties in predicting the behavior of complex systems 
is that engineers have trouble maintaining the 

9. Two economic aspects of a problem which most engineers are 
concerned with are and 



10 . Two types of information stored in the memory of a computer 
are and 



11. Modification of the by another set of operations 

is important for repetitive and logical abilities of a computer. 

12. The unit interrogates the memory for 
information as to what to do next. 

13. The gives the physical location of infor- 
mation in memory. 

14. An analog device computes by , a digital 

computer by applying a 
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15. The binary system is a number system of base 



16. Because the information in memory may be changed at any time, 
the digital computer is called a 

17. The physical property of is what makes a 

device an acceptable memory component. 
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SECTION E: THE FORTRAN LANGUAGE 



FORTRAN (FORmula TRANslation) is a language which allows the engineer 
to communicate a problem to a computer for solution. 

FORTRAN is not the natural language of a computer, nor is it the natural 
language of the engineer. Rather, it is a compromise between the two. 
To satisfy the computer, it uses symbols that the computer can understand 
and requires that the rules for their use be closely followed. To satisfy 
the engineer, it eliminates as many of the detailed computer control opera- 
tions as possible from the job of writing programs and uses a problem 
statement format close to that of the mathematical equation. 

The engineer describes his problem in the FORTRAN language; what he 
writes is translated into the natural machine language of the computer to 
be used in obtaining the solution. The translation is accomplished by the 
computer itself with the aid of a program called the FORTRAN Processor. 
The resulting machine-language program is then ready to be used to obtain 
the solution. This diagram illustrates the procedure: 



Problem 



COMPUTER 



... in FORTRAN 
language 




FORTRAN 
Processor 
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Problem 


... as a machine- 
language program 




\ 


t 










Results 




COMPUTER 


Input 
data 


W 


W 















Statements are the sentences of the FORTRAN language . They may: 

1. Define the arithmetic steps which are to be accomplished by the 
computer. 

2. Provide information for control of the computer during the execution 
of the program . 

3. Describe input and output operations necessary to bring in data and 
write the results. 



14 



4 . Specify certain additional facts such as the dimensions of a variable 
which appears in the program with subscripts . 

Because the printer and typewriter on a computer print only in upper-case 
letters, have a limited number of different characters, and are incapable 
of properly showing exponents and subscripts , FORTRAN statements at 
first appear somewhat confusing. 

An example of an arithmetic statement as it would appear on a FORTRAN 
coding sheet is: 

ROOT - ( -B+SQRTF(B**2-4.*A*C))/(2.*A) 

Translated, this says: 

The quantity to be known as root is equal to — that is , 
can be determined by — evaluating 



■VB2 



-B+VB'- 4AC 
2A 

where A, B, C are given values stored within the 
computer. 

Arithmetic statements look like simple statements of equality. The right 
side of all arithmetic statements is an expression which may involve 
parentheses, operation symbols, constants, variables, and functions, 
combined in accordance with a set of rules much like that of ordinary 
algebra. The symbols + and - are employed in the usual way for addition 
and subtraction. The symbol * is used for multiplication, and the symbol 
/ is used for division. The fifth basic operation, exponentiation, is rep- 
resented by the symbol **. A**B is used to represent A to the exponent 
B (that is, A B ). 

The FORTRAN arithmetic expression 

A**B*C + D**E/F - G 

v 
will be interpreted to mean 

A B C + D E - G 



That is, if parentheses are not used to specify the order of operations, 
the order is assumed to be: 

1 . exponentiation 

2. multiplication and division 

3. addition and subtraction 

Parentheses are employed in the usual way to specify order. For example 
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(A(B + C)) D 

is written in FORTRAN as 

(A*(B + C)) **D 

There are just three exceptions to the ordinary rules of mathematical 
notation. These are: 

1. In ordinary notation AB means A»B or A times B. However, AB 
never means A*B in FORTRAN. The multiplication symbol cannot 
be omitted. 

2. In ordinary usage, expressions like A/B-C and A/B/C are considered 
ambiguous. However, such expressions are allowed in FORTRAN and 
are interpreted as follows: 

A/B*C means (A/B)*C 
A*B/C means (A*B)/C 
A/B/C means (A/B)/C 

Thus, for example, A/B/C*D*E/F means ((((A/B)/C)*D)*E)/F. That is, 
the order of operations is simply taken from left to right, in the same way 
that 

A+B-C+D-E 

means 

(((A + B) - C) + D) - E 

3. The expression A B is often considered meaningful. However, the 
corresponding expression using FORTRAN notation, A**B**C, is not 
allowed in the FORTRAN language. It should be written as (A**B)**C 
if (A B ) C is meant, or as A**(B**C) if A( bC ) is meant. 

Besides the ability to indicate constants (like 3.57 and 2.), simple vari- 
ables (like A and ROOT) , and operations (like - and *) , it is also possible 
to use functions. In the previous example, SQRTF( ) indicates the square 
root of the expression in parentheses. 

Since the number of possible functions is very large, each computing cen- 
ter will have its own list of available functions , with information about 
their use. Functions given in this list must be referred to exactly as 
indicated . 

Some functions which might appear in a typical list are: 



16 



FORTRAN Symbol 

ABSF(X) | X | (absolute value of X) 

SQRTF(X) -v/x~ 

SINF(X) sinX 

COSF(X) cos X 

ATANF(X) arctan X 

EXPF(X) e x 

LOGEF(X) lo S e X 

Input output statements are those used to bring data into the computer to 
be stored for processing and to send results out. Typical examples are: 

READ 1, A,B, C This statement would cause the next card in 

the card reader to be read and the three 
numbers on it stored in locations assigned 
to the values A, B and C. 

PRINT 2, ROOT This statement would cause the number in 

storage identified as the variable ROOT to 
be printed. 

PUNCH 4,SUMA, SUMB This statement would cause the two values 

SUM A and SUM B to be punched on a single 
card. 

Similar input/output statements are included for reading and writing on 
magnetic tapes and magnetic drums and for such operations as rewinding 
or backspacing tapes. 

The numbers which follow the statements in the examples above specify 
the format in which the input or output should appear. Usually, a few 
standard formats will be specified for use in an installation. If a pro- 
grammer wishes to use something other than a standard format (for 
example, five numbers per card, in specified columns and in a specified 
notation), a new format can be introduced by means of a format statement. 
The FORTRAN Processor interprets the format statement and causes the 
machine -language program being generated to conform to the desired 
format in handling input or output statements which refer to it. 

Control statements allow the programmer to state the flow of his program 
through use of the logical powers of the computer. In writing a FORTRAN 
program , any statement in the program which is referred to by another 
statement must be given an identifying number (every statement can be 
given a number if desired but this is not necessary). The control state- 
ments refer to these identifying numbers for the purpose of branching 
from one part of the program to another. 

Important statements in this category are illustrated by the following 
examples: 
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GO TO 4 This statement indicates that the next 

statement to be executed (after having been 
converted to a machine-language program, 
of course) is the one numbered 4. 

GO TO (4, 18, 20, 40), K This statement is referred to as a computed 

GO TO since the value of K is computed in 
a previous statement. If, at the time this 
GO TO were executed, K = 3, then the third 
alternate (statement 20) would be the one 
chosen. 

IF (A**2 - B) 10, 20, 30 This statement allows one of three alternate 

next instructions to be chosen according to 
whether the expression in parentheses is 
less than, equal to, or greater than zero 
(in that order) . 

DO 25 I = 1, M, 2 This statement is a command to execute 

the sequence of statements which follow, 
up to and including 25 (or any specified 
statement number) , repeatedly. The first 
time, let the subscript I on the variables 
appearing in the statements equal 1 (the 
first number or variable following the = 
sign) . The second and each subsequent 
time the sequence is executed, increase 
I by 2 (the third number or variable following 
the equal sign — this can be left blank if it 
is to be 1) . When I becomes greater than 
M (the second number, or variable as in 
this case) do not repeat, but continue by 
executing the statement after 25. 

END This is the last statement in any FORTRAN 

program. 

In the DO example, reference was made to the changing of subscripts on 
variables. For example, A(I, J) would refer to a two-dimensional array 
of quantities a^. More will be said about this in the next section. 

What has been said about the FORTRAN language so far will be more 
meaningful by examining a complete FORTRAN program . The complete 
FORTRAN program to find the roots of the quadratic equation 

ax 2 + bx + c = 

might be written in the following logical steps: 

Arithmetic Statements 

The two roots are evaluated as arithmetic statements. A separate state- 
ment evaluating the discriminant is used to avoid the extra computation 
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that could be involved if this expression were written in both root state- 
ments and to permit later testing of the discriminant. 

2 D =B**2 - 4.0*A*C 

4 ROOT 1 = (-B+SQRTF (D) )/ (2.0*A) 

5 ROOT 2 = (-B-SQRTF (D) )/ (2.0*A) 

All the parentheses in statements 4 and 5 are necessary. D is enclosed 
in parentheses to define the expression to be operated on by the square 
root function. -B + SQRTF (D) is in parentheses so the additions and sub- 
tractions will be performed prior to the division. 2.0*A is in parentheses 
so the multiplication will be performed prior to the division. 

Input/Output Statements 

Statements are added to enter values for the constants A, B and C and to 
print the results of computation. Logically the read statement must pre- 
cede the arithmetic statements so that values will be available for compu- 
tation in the arithmetic expressions . The print statement must be after 
the arithmetic statements, so that values will be available for printing. 

1 READ 10, A, B, C 

2 D = B**2 - 4.0*A*C 

4 ROOT 1 = (-B + SQRTF (D) ) / (2.0*A) 

5 ROOT 2 = (-B - SQRTF (D) )/ (2. 0*A) 

6 PRINT 20, A, B, C, ROOT 1, ROOT 2 

In addition to the roots, the values of the constants A, B and C are printed 
to identify the results with the input. 

Control Statements 

Statements are now added for decision making and flow of the program: 
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Statement 3 is added to test for the mathematical possibility of a negative 
discriminant (complex roots). If the value of D is zero or positive, the 
flow of the program will be the normal flow to statement 4. If the value 
of D is negative, the flow will be to statement 8. Several alternatives are 
possible at this point depending upon the desired results: 

1. If complex roots are desired, they may be calculated and printed in 
their real and complex parts. 

2. This particular set of data could be ignored completely by transferring 
control to statement 1, where the next set of values would be read. 

3. Statement 8 could be a STOP statement which would cause the computer 
to stop execution of the program and permit manual intervention. This 
is not a recommended procedure . 

4. Statement 8 could be a PRINT statement to signal an error condition. 
In this program values A, B and C will be printed without the root 
values and can later be analyzed to determine what was wrong with 
this set of data. 

Statements 7 and 9 are added to complete the flow of the program by re- 
turning to statement 1 to read the next set of values. 

Statement 11 is the END statement to indicate the end of this program. 

This program will continue to be repeated as long as there are cards in 
the hopper of the card reader. When an attempt is made to read a card 
after the hopper has emptied, the computer will stop. 

Additional examples of FORTRAN programs will be found in the sections 
that follow. The use of additional statements will be explained at that 
time. 
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Exercises 



The FORTRAN language is designed to simplify engineering use of IBM 
computers. The use of the language with a particular IBM computer re- 
quires the use of a FORTRAN Processor for that machine. 

The power of IBM computers increases as the size and speed, in terms of 
machine components, increase. For the more powerful machines the 
FORTRAN language is expanded to take advantage of the additional types 
of components . For example , computers with magnetic tape input/output 
units allow use of FORTRAN statements to control the use of tape. 

Should the reader wish to become more expert on the FORTRAN language 
in general, or as it applies to a particular IBM machine, he may contact 
the local IBM representative for information on available IBM manuals 
and education courses . 



Write FORTRAN-language programs to: 

1. Put the numbers 1, 2, .0396, 30, 000.0 and -6.5 into the storage of 
a computer. 

*2. Compute C=A+B and print the resultant value for C, where A=3. 5 
and B=2.694. 

3. Compute C=A+B, where A and B are to be read into the computer and 
C is to be printed. 

*4. Compute C=A+B if B > 
C=A-B if B < 

where A and B are read into the computer, and C is to be printed. 

*5. Generate the integers 1 to 100 and print. 

6. Generate the value of sin x for x = .0000, .0001, .0002, 

.5. Print the value of x and sin x at each interval. 



*7. Change the program for solution of a quadratic equation to include 
solution for complex roots . 

* Answers to these exercises are found in Appendix I. 
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SECTION HI: SOME NEW AND OLD CONCEPTS 

Each new area of technology seems to construct a particular vocabulary 
for the necessary descriptions of efforts and communication of ideas. 
Computing is no exception to this rule; new terms such as "stored program" 
and "bit" have been added and new meaning given to old terms such as 
"iteration," "algorithm" and "memory." 

As far as use of computers is concerned, this growth in vocabulary is 
brought about by the new man-machine relationship. Control of a com- 
puting machine requires complete accuracy. This demands conciseness 
in statement of what is to be done, and preconception of what alternatives 
the computer takes as the result of all possible logical decisions it is 
asked to make. Of course, the usual give and take is allowed in regard to 
the testing of ideas . 

This section is designed to explain some of the procedures or techniques 
used in problem definition to simplify communication with the computer. 
In some cases this amounts to a re-emphasis of standard mathematical 
techniques . 

Block Diagramming 

The old cliche — a picture is worth a thousand words — has particular 
truth for the engineer in preparing a program for a computer . Since the 
engineer normally finds many uses for pictorial representations in his 
work, the block diagram is readily accepted as a powerful tool in comput- 
er programming. 

A block diagram is a picture of the steps which must be performed to ac- 
complish a particular job. The major function of the diagram is to clarify 
what must be done as a result of each decision. All the alternatives are 
accounted for in each case. In programming a computer it must be remem- 
bered that the computer has no way of anticipating the requirements of a 
problem and must be provided with all the information needed to reach a 
solution. 

The amount of information put into the block diagram, however, will de- 
pend upon personal preference, the programming techniques used, and 
the type of computer to be used. For example, Figure 1 illustrates two 
ways of block-diagramming the solution of a quadratic equation. The de- 
tailed diagram shows every step in the computation and is the type that 
might be used if the programmer were forced to work in machine language . 
The second diagram shows only the three main steps in the program but 
is quite sufficient for a programmer using the FORTRAN language . 

To further illustrate, consider a problem which involves the repetition of 
the same operation over and over. Let the problem be to create the sum 



Z = J^ pq - Yj 



This might be block-diagrammed in either of the two ways illustrated in 
Figure 2. 
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Read 
A, B, C 



AxC 



4xAxC 



A. 
DETAILED 




B*B 



B'' - 4AC 



D =\/B 2 _ 4AC 



D- B 



A + A 



X = 



D-B 

2A 



Punch X 



Figure 1. 
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Read X., Y. 
i=l,2,..,S 



ReadXj, Y- 
i=l,2, ..,5 



D 1 =[X 1 - Yj) 



D 2 = (X 2 -Y 2 ) 



D 5 = (X 5 -Y 5 ) 



Z = Dj + D 2 
+D 3 +D 4 +D 5 



Set i = 1, 
Z = 



Compute 
Z=Z + (Xi-Yj) 2 




© 



Punch Z 



Set i = i+1 



€) 



Punch Z 



Figure 2. 
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Subscripting 



The two diagrams illustrate an important programming concept — that of 
looping. Diagram A corresponds to a program in which each (X^ - Y.) 2 is 
computed separately and then all are summed to obtain Z . With only 
five values of i, this is not too unlikely a method of approaching the prob- 
lem. But, consider a situation in which i = 100 or 1000 or may vary 
according to some other characteristic of the problem of which this com- 
putation is a part. In such a case not only would the diagram be large and 
time-consuming to draw, but the program would also be large and time- 
consuming to write . 

In diagram B, advantage has been taken of the computer's ability to make 
logical decisions and to modify its own program . Since the number of 
times the computation is repeated depends only on the value of i (and the 
presence of successive values Xj and Y.), this one diagram — and the 
program which would be based on it — applies for any value of i. 

The algebraically incorrect statement i = i + 1 means replace the current 
value of i with i + 1. Similarly, Z = Z + (X^ - Yj) 2 means replace the par- 
tial sum, Zi_-p with the next value Zj. This use of subscripts brings us 
to our next topic. 



The purpose of subscripting in mathematics is to simplify notation and at 
the same time to remove ambiguity of meaning. In computer work the 
subscript is used in two ways: 

1. Space wise — to specify elements of arrays such as: 

A l A ll A 12 • ' ' • A 1N 

A 2 A 21 A 22 .... A 2N 

A 3 or A 31 A 32' ' ' * A 3N 



A M A M1 A M2 . . . A MN 

This allows reference to elements of an array through simple manipulation 
of the subscripts — just as in normal mathematical thought. 

2 . Timewise — to specify the chronological order in which a procedure 
occurs. For example, in Figure 2, diagram B, the subscript i is used 
not only to denote which member of the X and Y array is being operated 
on, but also to indicate how many times the computational step Z = Z + 
(X. - Y.) has been performed. This connotation and use is not so familiar 
to the engineer. The importance of the timewise aspect will become clear 
when iteration is discussed. 
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Absolute Value 



The absolute value becomes important because of the way in which a com- 
puter makes logical decisions. Computer decisions are made on the basis 
of the size of numbers — or, more specifically, upon whether a number is 
positive, zero or negative. For example, in order to do one of two things 
dependent upon whether A is (1) less than 500 or (2) equal to or greater 
than 500, the first thing to do is to subtract 500 from A: 

A - 500 = E 

Then the decision is based upon the size of E. Specifically, if E is neg- 
ative, procedure 1 is done, and if E is positive or zero, procedure 2 is 
done . 

If the difference represents the error in a procedure (that is, 500 is the 
true value and A is the estimate), the absolute value of the error must be 
less than a prescribed amount e and the test is: 

if I E | - e > do procedure 1 

or if |E | - e < do procedure 2 

The absolute value must be used, for in general there is no prior knowl- 
edge of whether the difference A-500 will be positive or negative. All 
the alternatives must be spelled out to the computer in its program. 

In many practical cases the "relative" error is a better measure of the 
error than the absolute error of a result. The relative error test for 
the above situation would be stated: 



if 

if 

Floating Point Arithmetic 



E 



500 
E 



500 



- e > do procedure 1 

- e < do procedure 2 



Engineering calculations often involve handling of very large numbers 
such as Young's modulus, which is 30,000,000 lbs. per sq. in. The 
common procedure is to write such a number as 3. x 10? to facilitate 
computation. 

A similar situation exists in handling very large or very small numbers 
in a computer. To get away from carrying a great many digits and to 
eliminate the effort of keeping track of the location of the decimal point 
for these numbers, a floating point notation is used. A common procedure 
is to maintain seven or eight most significant digits of a number plus a 
two-digit "characteristic" to indicate the proper position of the decimal 
point. The characteristic is developed from the exponent of 10 (assuming 
use of the decimal system) . 

For example, the number 

- .00000061957533 

26 



Errors 



can be represented as 

- .61957533xl0" 6 

However, this representation requires carrying two signs: one for the 
exponent and one for the fraction. This is inconvenient for computers which 
have only one sign associated with a storage location. Therefore, a com- 
mon practice is to add the exponent to some positive base . 

Using a base of 50, the characteristic becomes 50 + (-06) = 44 and the 
internal computer representation of the number becomes: 

4461957533- 

A FORTRAN programmer need not be concerned with this internal rep- 
resentation. The number might appear in a FORTRAN program as 

- .61957533E-6 



There are several sources of errors in computation which deserve con- 
sideration in any problem . 

Initial error: If x is the true value of a data reading and x* is the reading 
used in computation (reflecting an error in measurement, perhaps), the 
initial error is x - x*. 

Rounding error: This type of error results when the less significant digits 
of a quantity are deleted and a rule of correction is applied to the remain- 
ing part. For example, pi, 3. 14159265. . . , rounded to four decimals, 
is 3.1416. 

Truncation error: To simply chop off at four decimal places for pi, giving 
3.1415, would result in a truncation error. Another common source of 
truncation error is in chopping off all terms in an infinite series expansion 
after a particular term. For example, cutting the series for e x at 

e x = 1 + x + x 2 + x 3 
"21 3! 

gives truncation error — sometimes called residual error for series 
approximations . 

Propagated error: If x is the true value of a variable and x* the value 
used in computation, then f(x) - f(x*) is the propagated error. Errors 
may build up in computation! 



Recursion Formulas 



Very often computation of a set of numbers can be simplified by formulas 
which give one member of the set in terms of other members of the set. 
An example of this is finding the new coefficients of a polynomial when the 
original polynomial has been reduced by the extraction of a root. If r is 
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Algorithms 



one root of the polynomial equation a n x n + a ^ x 11 " 1 + a^_o x 11-2 + . . . 
a = 0, then the origina 
the reduced polynomial: 



a = , then the original polynomial can be divided by (x - r) to obtain 



V*" 1 + (a^ + ra n ) x^-2 + ^ n _ 2 + r <a n _ 1 + r^)] x*"3 +. . . . 



c-r^T?- - - 11 - 1 - - " n " 2 



+ a x" x + a x 11 " +. . . 
n n-1 n-2 



a x n - ra x 11 " 1 
n n 



<Vi + r V xn_1 + *-2 xn " 2 + * ' ' ' 

(a n _ 1 + ra n ) x 11 " 1 - r (a n _ 1 + ra n ) x 11 " 2 



[ a n-2 + r < a n-l + r V] x 



n-2 + 



k k-1 
Let: Cfcx + c k _jx + + c-jx + c Q = 

be the resulting polynomial (where k = n-1). 
Then, equating coefficients gives: 



Ci = a 
k n 

c k-l = < a n-l + ra a) = an-1 + rc k 

c k-2 = [ a n-2 + r < a n-l + ra n>] = a n-2 + rc k-l 



or in general 



c. = a. . + re. „ 
i l+l l+l 



This is a recursion formula for c. . This formula has obvious advantages 
for programming the calculation of the new coefficients . 



An algorithm is a theorem which may state that a solution to a problem 
exists, and/or a procedure for obtaining the solution. The term is fre- 
quently encountered in computer literature because the form of an equa- 
tion is often all-important in programming efficiently. 

Two examples will be given here: 

1. To compute the value of the polynomial 

f (x) = a ft x 5 + a-j x^* + a^x^ + a,x 2 + a. x + a 5 
computation is simplified if it is written as 

f(x) = (((((a ) x + a x ) x + a 2 ) x + a 3 ) x + a 4 ) x + a 5 

which indicates that the polynomial may be computed by repeating 
the steps "multiply and add" five times. 
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2. The square root of a positive real number, A, can be computed by 
using the algorithm 



1+1 =T< x i + |:> 



This example also introduces the idea of iteration. What is meant 
by this equation is that to obtain an estimate of A , start with a first 
guess x ft . Substitute it in the right side of the equation to obtain the 
next estimate of A* . A few of the steps are: 

x i=| (x o + # > 



= ±< v + A 
: 2 



x 3 = i( x 2 + £> 



x. , -. = — (x. + — ) 
i+l 2 v i x. 
l 

Note that the algorithm states a procedure for solution of the prob- 
lem, not just one formula evaluation. 

Further, note the timewise use of the subscript i; i + 1 indicates a 
result dependent upon the previous result subscripted by i. 

Taking a value of A (say 25) for which the square root is known, and 
performing the indicated operations will make this algorithm clear. 
If 2 is used as a starting estimate, then; 

x o = 2 

X l = 2~( 2+ ^ =7 ' 25 



1 / 2b \ ^ 

x =—(7.25 + ) =5.35 

2 2 V 7.25/ 

X 3=T( 5 - 35+ 5f 5 )^- 01 



There are several ways of deciding the point at which to discontinue 
the iteration: 

a. When differences between successive estimates are less than 
some £ of accuracy | x. + , -x. | < £ 

b. When relative differences between successive estimates are 

less than some £ of accuracy x^ +1 -Xj < £ 

x. 
l 
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This makes the choice of £ independent of the size of the original 
function value A. 

c. When relative differences between the function of the estimate and 
the original function value A are less than some £ of accuracy 



(x i+1 ) -A 



This eliminates the possibility of an erroneous solution in a sit- 
uation where convergence is slow and is the preferred type of 

test. 

The general statement for this test in finding a value of x which 
satisfies 

f (x) = A 



is 



f (x i+1 ) - A 



In succeeding chapters the problem of solving for values of x which 
satisfy equations of the form x = f (x) are discussed. Methods a and 
b remain the same but the test equation corresponding to c above is 



x. 



i+1 



" f < x i+l> 



f (x i+1 ) 



Complex Number Calculations 



Complex numbers can create some questions of handling and yet are of 
great importance to electrical and electronic engineers . The following 
rules handle complex arithmetic: 

A. Addition 

(a+bi) + (c+di) = (a+c) + (b+d) i 

B. Multiplication 

(a+bi) (c+di) = (ac-bd) + (ad+bc) i 

C . Division 

a+bi = (a+bi) (c-di) = (ac+bd) + (bc-ad) j 
c+di (c+di) (c-di) (c 2 +d 2 ) (c 2 +d 2 ) 
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D. Square Root 

If we let u+vi = (a+bi) 2 

squaring both sides gives u 2 + 2uvi - v 2 = a+bi 

Equating the real and imaginary parts 

(1) u 2 - v 2 = a 

(2) 2uv = b or v = i 

v ' 2u 



then substituting from (2) into (1) 

b 2 , 
4u2 



u 2 - b 2 - a = 



± 

' , 2 



2 
and multiplying through by u gives 

u 4 - au 2 - b 2 = 
and solving the quadratic gives 

Only the positive sign is considered since u is real; thus 

mu -(* + (@' + (*)')** 

If a is positive, compute u by (3) and use (2) to find v. 

If a is negative, compute v by (4) and use (2) to find u. 

E . For a general power of a complex number use 

(a+bi) n = rV 11 " 9 " where r = (a 2 +b 2 ) T 

■e=tan~ 1 A 
a 

and e = (cos -9- + i sin-6-) = cos n-6-+ i sin n0 

then (a+bi) n = r n (cos n-9- + i sin n-6- ) 

Transformations and Identities 

Recognition of mathematical transformations and identities which are 
possible may simplify computation. 
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For example, assume in a particular problem that the value for the tan- 
gent of an angle is given: 

tan •&■ = A 

and that cos -9- is to be computed. One approach would be to find 
-0- = tan" ■'■A and then compute the cos^- . However, the cosine may 
immediately be computed as 



cos -9- = 1 



\ll + A 2 




Exercises 



saving the evaluation of two trigonometric functions . 

A useful transformation in computing is the logarithmic transformation. 
Suppose we have an exponential function y = k x a . If we take the natural 
log of both sides, we have 

l°g y = a l°g x + l°g k 

thus reducing an exponential relationship between x and y to a linear rela- 
tionship between the transformed variables log x and log y. 

Most of the mathematical techniques which have been outlined here are 
quite familiar to the average engineer. However, they are given this 
emphasis because recognition of their use in particular situations can 
greatly simplify solution of computer problems . 



1. For the set of tabular data shown, construct the difference table to 
determine to what degree of x is y dependent. 









1 


3 


2 


18 


3 


57 


4 


132 


5 


255 


6 


438 



2. Divide 8 + 12i by 3 + 2i. 



3. Using the iteration formula for finding the square root of a number, 
find the square root of 12 accurate to 7 places (the absolute difference 
between two successive estimates of 412 is less than .00000005). 

*4 . Using only the functions SINF (X) , C(j)SF (X) , ATANF (X) , EXPF (X) , 

L(J)GF(X) and SQRTF(X), write the FORTRAN statement for computing 
all trigonometric, inverse trigonometric and hyperbolic functions. 
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5. Construct a block diagram chart to show the logic of computing the 
value of the polynomial 

f (x) = a 5 x 5 + a^x + a 3 x + a^x 2 + a-jx + a Q 

where the values of the coefficients a., , a , a are to be 

1 2 6 

read into the computer and the value of f (x) is to be punched out. 

*6. Repeat 5 and generalize the block diagram to show the logic of han- 
dling any degree polynomial. That is, block-diagram the computation 
of f (x) , where 

f (x) = a-|X n + a 2 x n 4- a n+1 

7. Write a FORTRAN program for exercise 6. 

8. The block diagram in Figure 2 is incomplete. Why? 
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SECTION IV: PRINCIPLES OF ITERATION 



Arithmetic, Mathematics and the Computer 

Although arithmetic is one class of mathematical operation, for purposes 
of clarity it is best to differentiate arithmetic operations from the other 
mathematical operations in discussing computers. Arithmetic is performed 
when any of the operations of addition, subtraction, multiplication and 
division are applied to a pair of numbers . Other mathematical operations 
are performed when one of a variety of operations such as algebraic sub- 
stitution, integration, differentiation, expansion to series representation, 
etc. , is applied to one or more symbols. In essence arithmetic deals 
with numbers while mathematics is the manipulation of symbols. 

Computers are designed to efficiently perform all the arithmetic opera- 
tions. Most computers can also handle in a limited way the alphabetic 
and some special characters which might be used as symbols (as is obvious 
from the use of the computer to translate FORTRAN programs). 

The computer, as designed, cannot do original mathematics. However, 
once the mathematics of a problem are accomplished, the computer is 
most useful in obtaining the required numerical values by performing the 
necessary arithmetic operations. 

In addition, once a particular mathematical problem has been properly 
solved and the program for arithmetic evaluation on a computer has been 
completed, all similar problems, differing only in data values, can be 
henceforth evaluated by the one program. Two examples will illustrate 
these facts: 

A. Given the two simultaneous equations 

a-i -jX-i + a-i <y&-o ~ "i (-*-) 

a 21 X l + a 22 X 2 = b 2 (2) 

find values of x and x which satisfy these equations when 

a n =2, a 12 =2, a 21 =3, a 22 =4, b^S, b 2 =6 

To solve this problem we have distinct mathematical and arithmetic steps . 

Mathematically: 

Algebraically using (1) write x, as a function of x 2 

x =f fc ) = b l " a 12 X 2 
11 

Substitute this function f-, (x,) in the second equation to obtain an equa- 
tion in Xr, 
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_a 21 a 12 x 2 + a 21 b l + a 22 x 2 = b 2 



(4) 



u ll 



u ll 



This equation may then be treated algebraically to give x„asa function 
of the coefficients only 



x 2 = f b 2 - a 21 bj \ / / a 22 - a 21 a 12 



)/( 



11 



11 



Arithmetically: 

Substitute the values given for the coefficients in equation (5) to deter- 
mine x 



x = /6-3.5 
2 T 



t-3 _i_ 2 ^ =-1.5 
2 



Use this result in equation (3) to find a value for x 



x a 5-2(-1.5) =4 
1 2 



The computer solution to this problem considers only the arithmetic op- 
erations implied by equations (3) and (5). A FORTRAN program to eval- 
uate these two equations is shown in Figure 3. 



I STATEMENT 
t NUMBER 

1 5 


i 
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Figure 3. 

NOTE the following items of interest concerning this program: 

1. Only the solution equations have been programmed. The original 
equations are not of interest to the computer since the mathematics 
must be performed before using the computer. 

2. Equations (3) and (5) are programmed in their symbolic form and actual 
values for the coefficients are read into the computer before compu- 
tation takes place. By programming in this way the computer can now 
evaluate the solution of any set of two linear equations in two variables, 
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in addition to the particular example discussed. In a sense, by our 
programming, we give the computer the capability of doing nonoriginal 
mathematics. 

B. Solve the following differential equation with stated initial conditions 
for y when x is equal to 1: 



dy 

dx" = x y 



(1) 



y = 1 when x = 
Mathematically : 

Integrate (1) for y as a function of x 
dy 

dx- = x y 



dy 



x dx 



lny =— +k 

y = e * 

Considering the initial condition y = 1 when x = 0, then 

k = 

Arithmetically: 

x£ 

Since the analytic solution is y = e 2 it is now necessary only to sub- 
stitute the value of x for which y is to be determined and perform the 
indicated arithmetic . 

y =e 2 = 1.64872 

A FORTRAN program to accomplish the required calculations is shown 
in Figure 4. 
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Figure 4. 
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NOTE the following: 

1 . Again only the analytic solution is programmed . 

2. The program is written so as to allow many values of y to be com- 
puted over any range in x . 



Numerical Analysis 



The two problems chosen happen to be of types for which the mathematical 
analysis is well defined. The variables in question may be equated to 
functions of the known quantities — analogous to the formula y = f (x, z) as 
discussed in the first chapter. 

Various analytic techniques successfully handle a large class of the math- 
ematical problems which arise in the description of the physical world. 
However, a vastly larger number of mathematical problems are not sub- 
ject to the normal analytic procedures of algebra, analytic geometry, or 
the calculus. 

For the theoretical mathematician or physicist it may be sufficient, when 
meeting these problems, to be able to show that one or more approxima- 
tions will allow an analytic solution of the problem to any desired degree 
of accuracy. 

The engineer, on the other hand, must have a solution from which accurate 
numerical results can be obtained if the theory is to be of value. These 
results must be obtained at a reasonable expense of time and effort. A 
mathematician may be content to know that a series will converge. An 
engineer must be concerned with the number of terms in the series re- 
quired for a given accuracy. 

Numerical analysis provides techniques for finding concrete numerical 
results for mathematical problems of the type under discussion. This 
manual describes many physical problems which require numerical analy- 
sis. Some of the simpler analysis techniques are discussed. The reader 
is referred to the many excellent texts available on numerical analysis for 
an adequate discussion of particular analysis theory. 



Algebraic Equation 



Problem A belongs to a very special class of algebraic problems, a set 
of n nonhomogeneous linear equations in n variables. Fortunately, a 
large class of engineering problems such as the force and moment prob- 
lems of stress analysis give rise to this simply handled type of equation. 
The diagram below shows a simple force equation example giving rise to 
a set of two linear equations . 
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Zr y = o 




T cos-0-- T cos 0=0 

T 1 sin-9-+ T 2 sin (J) = W 

T-. and T are readily solved for specified values of -9- , (|) and W. 

Consider a set of equations similar to A but containing a simple 
nonlinearity: 



a ll x l +a 12 x 2 ~ b l 

annXn 

a 21 X l + G -\ 



(1) 
(2) 



The previous algebraic steps would give: 



a _12_. , bi 
a ll " a n 



x„ = - — x 2 + ^_ (3) 



and: 



a ll 11 

or with appropriate substitutions for (4) 
bx 2 



x 2 = ae 



+ c 



(5) 



(4) 



This gives an equation in x 2 only, but from the equation it is impossible 
to directly evaluate x 2 - One approach is to expand e bx2 to several terms 
of its series: 

2 3 
e bx 2 = ! + bx 2 + (bx 2 ) + (bx 2 ) + 



2! 



3! 



However, for sufficient accuracy it may be necessary to carry a large 
number of terms in the series resulting in a polynomial equation being of 
degree greater than 4. In such a case, the equation is again impossible 
to evaluate directly. This chapter will discuss handling equation (5) and 
polynomials by means of iteration. 
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Equations of Calculus 



Similarly the differential equation of example B has an analytic solution. 
However, the differential equation 



dy 
dx 



= * + y 



has no corresponding analytic solution. To find values of y for corre- 
sponding values of x, a numerical approach is necessary. Likewise, 
many integrals may be obtained only by numerical approximations. Again, 
some simple techniques will be given for handling this type of problem. 



Iteration Type Problems 



There are three principal reasons for the use of iteration with computer 
handling of engineering problems . 

1. The mathematical statement of the physical problem requires an itera- 
tion approach for evaluating one or more of the variables . 



The chart below shows examples of four types of mathematical problem statements. 



Single Equations 


Multiple Equations 




i. x/2 

x= 5 e 


M H =,A H -m 

M He = A He- 2m 

R„ - R 
He H 


Direct Iteration 
Sufficient 




R H/M H - R He/M He 




= tan0-K 


B = 90°-*-sta" 1 ( L + * + E ) 

L = A sin (90° - B - % -<£) - (R-E) 


Numerical Techniques 
Required 



2. Many designs are to be evaluated in a search for the best design, that 
is, optimization of design. Often this problem can be reduced to the 
previous problem by selection of an appropriately expressed mathemat- 
ical criterion of what is optional. 

3. The mathematical expression for the physical problem is to be eval- 
uated for many values of one or more known parameters — such as 
time, in a problem of motion, or degree of rotation, in a geometry 
problem . 



A. One problem which illustrates iteration is the solution of an equation 
which arises frequently in absorption problems of optics, electric 
fields and nuclear engineering: 



x = ae bx 



(1) 
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where a and b are constants . 

This cannot be solved explicitly for x, so the following iteration pro- 
cedure is used: 

1. Make a guess at x. 

2. Use this in right-hand side of equation (1) to give a new value for x. 

3. Call this new value of x the next guess. 

4. Repeat steps 2 and 3 until two successive guesses either agree or 
differ by an amount less than the allowable error. 

A FORTRAN program to solve such a problem, where 

x -- 0.2e 0,5x 

is as follows: 

READ 10, A, B, X 



1 PRINT, X 

2 XNEW=0.2*EXPF(0.5*X) 

TEST=ABSF(X-XNEW) 
IF(TEST-.00005)4,4,3 

3 X=XNEW 
GO TO 1 

4 PRINT, XNEW 
END 



Read values for A, B and 
initial estimate of x (1.0). 

Print estimate of x. 

FORTRAN arithmetic state- 
ment of equation to be solved. 

Find absolute value of dif- 
ference in last two estimates. 

If difference is less than or 
equal to .00005, go to step 
4. Otherwise, go to step 3. 

Store new estimate of x. 

Return to step 1 and repeat. 

Print last estimate of x. 

End of program. 



If this were translated into machine language and the resulting program 
run on a computer, the successive estimates of x would be: 

x 

1.000000 
.329744 
. 235848 
.225032 
.223818 
.223682 
. 223667 
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The last value in the list satisfies the requirement that the difference 
between it and the previous estimate be less than 0.00005 in absolute 
value. Hence it is the solution. 

B. The simple iteration procedure in problem A seems like a nice way 
to solve problems — and in many cases it is. However, on applica- 
tion to problem B — much like problem A except that a constant has 
been changed — it fails: 

f(x) = x - Ae Bx = A = .2, b=3 (1) 

The reason for the failure of convergence of the direct iteration technique 
should be explained. Direct iteration has the formula 

x _,_ = f (x ) 

n+1 N n 

This indicated procedure will converge only for 

-l<f (x)< 1 

Bx 

where f (x) = A e for this problem, and the prime indicates the first 

derivative. The choice of a larger B makes f ' (x) > 1 and therefore makes 
convergence by the direct procedure impossible. 

Since simple iteration fails to produce convergence, a different method to 
obtain an estimate must be found. The Newton-Raphson method answers 
this need. In this instance, successive estimates are obtained from 



F < x m> 
c m+l x m 



F ' < x m> 

From equation (1): 

Bx™ 
F ( x m) =x m- A e 

and 

F'(x m ) = l-BAe BXm 
Applying the Newton-Raphson formula: 



x - Ae 
= x - 



Bx m 



m+1 m 1 - BAe Bx m 

To solve this problem in FORTRAN, replace statement 2 in the previous 
FORTRAN program with the arithmetic statement for (3): 

2 XNEW=X- (X-A*EXPF(B*X)) /(1-A*B*EXPF(B*X)) 

If this FORTRAN program were translated into a machine-language pro- 
gram and run on a computer, the successive values of x would be: 
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1.00000 
. 19742 
.22364 
.22366 

Note the increased rate of convergence obtained using this new method; 
only three new estimates are computed, compared with seven in the pre- 
vious example. 

C . The following engineering problem will illustrate one origin of the 
type of equation just discussed. The analysis makes use of the clas- 
sical optimization principle of equating the derivative of a function of 
one variable to zero. The values of the variable which satisfy the 
resulting equation are those for which the original function is either 
a maximum or a minimum . (Since engineers are often interested in 
optimization of design, it is unfortunate that so few design problems 
can be expressed in a form suitable for the use of this principle.) 

For high-potential conducting through walls, tubular insulators are used 
which are covered with metal on the inner and outer surfaces. 



Bore 




> r <r 



What must be the ratio of the external diameter, 2R, to the bore width, 
2r, for the cross section Q to be a minimum? 

The ratio q of the line voltage to the maximum admissible field strength 
is given by 



q = r In R/r 

The cross section Q is given by 

9 9 
Q= tt-(R- - r ) 

According to (1), letting x = R/r, then 

r = q/ln x 



(1) 



(2) 



or 



Q = 7T q 2 (x 2 - l)/(ln x) 2 



(3) 
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Since Q is to be a minimum, it is necessary to find a point on the curve 
defined by (3) where the slope is zero. To do this we take the first deriv- 
ative and set it equal to zero: 



dQ 
dx 



■n-q 



2x 



2(x^ - 1 ) 

(lnx) 2 x(lnx) 3 



(4) 



From the bracketed expression: 



x - (x - l)/x In x = 



or 



In x = 1 - 1/V 
This in turn gives 



x - e 



1-1/x' 







(6) 



The solution for x is readily recognizable as the type of problem just 
discussed. 

D. Another example can be taken from an equation which occurs when 
working with helical gears . This is 



K = involute of (J) = tan (|) - d) 



(1) 



where K is found in a table and it is necessary to find an angle (J) satis- 
fying the equation. 

The Newton-Raphson method can again be used but it is worthwhile to stop 
and consider a good way to get a first estimate of (|) . Since the rate of 
convergence — and in some cases, the possibility of converging — will 
depend on the initial guess, this is a subject worthy of some thought. 

The following graphic presentation of this problem will clarify the situa- 
tion. 

Plotting Newton's function: F(0) = tan (J) - (|) - K gives the curve as shown: 




Possible range of 

•< — initial estimate — ^ 

for convergence 
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Enlarging the area of intersection of F ((|)) with the (|) axis shows the basis 
for Newton's procedure. 




If (J) is the first estimate of the root, the figure indicates that the next 
good estimate would be the intersection of the projected slope of F ( $ ) 
evaluated for ty n . The following use of the construction angle ■©• shows 
this to be precisely Newton's procedure: 



tan-9- 



F(d> n ) 



n+1 



= F' 



,) 



or 



d> = _F((D ) 
n+1 n P^TFl 



The dotted lines in the figure show the limits on the range of the initial 
estimate for which Newton's technique will converge to the proper value, 



A good first estimate can often be determined by a careful examination of 
the problem to be solved. In this case the first two terms of the trigono- 
metric series for tan (j) give: 
,h3 



tan (J) = (|) + 



3 



(2) 



Substituting this in (1) and solving for (J) gives 



<D ~ (3K) 



1/3 



Thus an initial guess which should be reasonably close to the solution can 
be computed from the given value s of K. This can Toe ah imporXanf advan- 
tage where convergence is likely to be slow. 

To illustrate the effect, the problem has been solved for two values of 
K (0.001 and 0.01). In each case, the initial estimate of (J) = 0.52359 
radians (30°) was used rather than a value computed as described above . 
The results obtained are as follows: 
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(K = .001) fl) (K = .01) 

.52359 .52359 

.36535 .39235 

.25482 .32546 

.18624 .30818 

.15296 .30702 

.14440 .30701 
. 14385 
. 14372 

Note that for K = .01 fewer iterations were required because the initial 
estimate was closer to the correct result. Had the approximation formula 
developed above been used, the initial estimates would have been 0. 14423 
and 0.31072 respectively, and no more than two iterations would have been 
required in either case. 

E . A simple geometry problem which lends itself to iterative solution is 
that of finding the angle for which the arc and chord in a circle of 
given radius will enclose a stated area. 




Area A is given by 

2 9 

A = 7r r -6- _ r* sin -a (i) 

360 2 

Since -0- cannot be solved for directly, the Newton-Raphson technique is 
used to find a solution. Converting -©- to radians and applying the N-R 
formula gives the iteration equation: 

-0- --a- _ F(-0- ) 
n+1 n _, .n . 

where 2 2 

F(Q)=- r » n - r sin ^ n " A 
2 2 

and „ 9 

F'(-e-) =£_ - r z cos-9- n 
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It is interesting to note how frequently an iterative problem is given rise 
to by asking an unfamiliar question to a familiar relationship. Generally 
equation (1) is used to determine the area A. In a sense, iteration arises 
when the "answer" is given and the "question" is to be determined. 

F. Very often iteration can be the most convenient procedure in finding 
a root of a cubic such as 



3,2 
ax + bx + ex + d 







Newton's technique will generally suffice here. 

For polynomials of degree greater than 4, some type of iterative proce- 
dure must be used. There are many techniques available to handle the 
higher-degree polynomials, but they will not be discussed here. One of 
the important advantages of the FORTRAN language is that it makes it 
more readily possible for computing organizations to establish standard 
programs for such problems and exchange them. 

G. This problem is an example of how a simple circuit can give rise to a 
cubic equation. Generally it is simpler to make successive approx- 
imations than to solve the cubic . 



Nonlinear 

Passive 

Component 

i^kv 3 / 2 

k = 10" amps/ volts 3 / 2 



T 



Linear 
Active 
Component 

Emf = 100 volts 

R = 5000 ohms 



Problem: Given the circuit shown, find the potential difference V from 
(1) to (2) and the current I flowing in the circuit. 

Using the basic relationship for the potential drop, 

V = E -IR 



and considering the active component , gives 

.3/2 



or 



V = 100 - 5000kV 
,3/2 



V = 100 - . 05V 
a cubic equation easily solved using Newton's technique. 
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H. An interesting example of the use of iteration can be taken from the 
computations employed in spectroscopy. Spectroscopy and astronomy 
are two areas of science in which it is meaningful to speak of numbers 
accurate to seven or eight significant figures . The success of iterative 
approximations in the field of spectroscopy demonstrates beyond ques- 
tion its possibilities for accuracy. 

Given experimentally determined values of the Rydberg constants for 
hydrogen (R H ) and helium (Rjj e ) and the atomic mass (A H , A R ) of 
the elements, it is possible to compute the mass of the electron. 

Assume the following data is given 

R H = 109,677.68 cm -1 A R = 1.00812 

R~ = 109, 722. 34 cm -1 A TT - 4. 00388 

He He 

if R^ is the Rydberg constant assuming the mass of the nucleus to be 
infinite, the following relationships exist where m is the mass of the 
electron: 

R co = < 1+ M R ) R H (D 

R = (1+— ) R 
oo * M He ' He (2) 

The atomic mass includes mass of the electron: 



Ajj = M R + m (3) 

A u = M tt + 2m 
He He 

This gives a useful first approximation 



M H = A R (4) 

He He 
Equating (1) and (2) gives a first evaluation of m: 
R„ R, 



m 



= He - H (5) 

R R 
H He 



M M 
H He 



To obtain subsequent approximations, equation (3) is rewritten as: 

M H = A R - m (6a) 

and 

M He =A He- 2m < 6b > 

and, using the new values m is re-evaluated by (5). 
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A FORTRAN program for this procedure is: 



20 EMASS=0.0 

30 RH=109677.68 

40 RHE=109722. 34 

50 AH=1. 00812 

60 AHE=4. 00388 

70 EPSLN=. 00000005 

1 EMH=AH-EMASS 

80 EMHE=AHE-2.0*EMASS 

90 EMASP= (RHE-RH)/(RH/EMH-RHE/EMHE) 

100 PUNCH 20,EMASP,EMH,EMHE 

110 IF(ABSF(EMASS-EMASP)-EPSLN)3,3,2 

2 EMASS=EMASP 
120 GO TO 1 

3 END 



Running this program on a computer produces the following results: 
_m_ M H M Re 



.00054871426 1.0081200 4.0038800 

.00054836568 1.0075713 4.0027826 

.00054836592 1.0075716 4.0027833 

It is significant to note that the simplest iterative technique produces the 
necessary accuracy for this application in only three iterations. 

I. Earlier it was stated that iteration could be used in obtaining numerical 
solutions to differential equations. To illustrate, consider a problem 
for which the exact integral is known: 

dy (i) 

dx y 

and find y as x varies from to 1.0. The initial conditions are that for 
x = 0, y = 1.0. 

Since discrete values of dx and dy must be used, it is necessary that these 
be replaced with the approximations: 

Ax~dx 
and 

Ay~dy 

If a fixed value is assigned to A x in the range to 1.0, then: 



X =0 



x l = x + Ax 
x 2 = x l + Ax 



x n+l =x n +Ax 



(2) 
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To find the values of Ay, equation (1) is changed to the approximation: 

Ay =x ..y Ax (3) 

^n n+l^n ^ ' 

Then, 



y „ = y + Ay 

J n+1 'n J n 



(4) 



Repeated evaluation of equations (2) , (3) and (4) will give the solution: 



3 


DELX= .05 


4 


X=0.0 


5 


Y=l .0 


1 


PRINT 20,X,Y 


6 


X=X+DELX 


7 


DELY=X*Y*DELX 


8 


Y=Y+DELY 


9 


IFCX-1 .0)1,1,2 


2 


END 



Taking A x = .05, the following solution was obtained 



x 



.00000 
.05000 
.10000 
.15000 
.20000 
.25000 
.30000 
.35000 
.40000 
.45000 
.50000 
.55000 
.60000 
.65000 
.70000 
.75000 
.80000 
.85000 
. 90000 
.95000 
1.00000 



1.00000 
1.00250 
1.00751 
1.01507 
1.02522 
1.03803 
1.05361 
1.07204 
1.09348 
1.11809 
1.14604 
1.17756 
1.21288 
1.25230 
1.29613 
1.34474 
1.39853 
1.45796 
1.52357 
1.59594 
1.67574 



The analytic solution to (1) for the stated boundary conditions is 

2 
x 

and for x = 1.0, y = 1.64872. The error in the numerical solution is 
1.67574 - 1.64872 = 0.02702, or less than 2%. Considering the size of 
the interval, this is not out of line. Increased accuracy could be obtained 
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by decreasing Ax or by using a more sophisticated technique — and there 
are many. 

Consideration of the first alternative should take into account that de- 
creasing the value of A x may increase the rounding error and hence make 
it necessary to carry along a greater number of significant digits . The 
second alternative must be viewed in the light of the accuracy required. 
The more sophisticated techniques will require greater programming 
effort and more computer time . It is not reasonable to spend this time 
and effort to obtain accuracy of .01% if the required accuracy is only 2% — 
particularly in cases where the input data is accurate to, say, 5%. These 
points should be considered in programming any problem for a computer 
if the most economical use is to be realized. 

J. It is easy to extend the simple numerical integration procedure just 
discussed to higher-order differential equations. For example, sup- 
pose the following is to be solved over the interval to 1.0: 

i?y + A-^ + Bxy = 
dx^ dx 

where x = 0, y = 1.0, (dy/dx) = 1, A = 2 and B = 3. 

The analytic solution must be in terms of an infinite series . Therefore , 
whether a series is determined or numerical integration is performed, 
the solution involves considerable computation. The series representa- 
tion allows control of the error while the stepwise integration procedure 
to be described may not. 

An approach to this problem is to start with 

^2J0 = " A(dy/dx) - Bx y 
and then obtain successive values of each of the variables from 

y n+l = y n + < dy/dx >n Ax 

X n + l =X n +Ax 

(dy/dx) n+1 = (dy/dx) n + (d 2 y/dx 2 ) n Ax 

< d Vdx\ +1 = -A(dy/dx) n+1 -Bx n+1 y n+1 
A computer solution to this problem is shown in Figure 5. 



* * * * 



In the first few problems just presented, iteration was used to improve 
the estimate of a single solution. In the case of the differential equations, 
iteration provided successive values of y for finite changes in x. These 
are two entirely different concepts of iteration. 
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1 A=2.0 

2 B=3.0 

3 XN=0.0 

4 YN=1 .0 

5 DYDX=1 .0 

6 D2YDX=-A*DYDX-B*XN*YN 

7 PRINT 20, XN,YN,DYDX,D2YDX 

8 IF(XN-1 .0)9,14,14 

9 YN = YN+DYDX-*.05 

10 XN=XN+.05 

11 DYDX=DYDX+D2YDX*.05 

12 D2YDX=-A*DYDX-B*XN*YN 

13 GOTO 7 

14 END 



XN 



YN 



DYDX 



D2YDX 



Exercises 
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Figure 5 

Engineering computational problems may be classified in the following 
way — the first category being of trivial interest here, although certainly 
not trivial as far as computing is concerned: 

1. Parameters specified, explicit equations . 

2. Parameters specified, iteration required. 

3. Parameters variable, explicit equations. 

4. Parameters variable, iteration required. 

Problems A through H are typical of the second category. Differential 
equations may fit in either the third or fourth category — examples I and 
J falling in category 3. The next subject to be discussed — mathematical 
models — also comes under either category 3 or 4. 



1. For problem C, find Newton's iterative equation and write the FORTRAN 
program to solve the problem. 
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2. Write Newton's equation to find a root of 



5x 3 + 2x 2 + 3x + 1 = 



3. The following is a problem which occurs in the planetary gear systems 
currently used in automatic transmissions on automobiles. The curves 
y = sin ait and y = e~ a * are plotted on the same set of y, t axes. De- 
termine the smallest positive value of t at which they intersect. 




a. Set up the equation for which a root is to be determined. 

b. Write Newton's iterative equation to find the root. 

*4. A common numerical procedure for calculating integrals is the use of 
Simpson's rule which may be stated: Given y as some function of x, 
the integral of the function y over the range x = x , to x = x may be 
approximated by: 



x = x„ 



'x = Xi 



~Ax 



y dx =^- (y x + 4 y2 + 2 y3 + 4 y4 + 2 yn _ 2 + 4^^ + y n ) 



where x^ = x^ +Ax 

and y. is evaluated for the corresponding x. . 

Write a FORTRAN-language program to evaluate 



S 



x = 1.0 



y^dx 



x = 
using a A x = .05. 



where y^ lT% 2 
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SECTION V: MATHEMATICAL MODELS 



It is just as well not to belabor the question of what makes a mathematical 
model. However, it might be defined as a set of relationships which define 
a physical system. Different sets of operating conditions or input para- 
meters, when handled by the model, give correspondingly different sets 
of operating results — problem category 3 . If iteration is required for a 
set of input parameters, the problem is in category 4. 



Solenoid Design 



This problem is most typical of design problems where several variables 
are involved and there is no obvious procedure for arriving at the best 
design (category 4 — parameters variable, iteration required). 

Problem Statement: 

Choose a wire for a solenoid such that the power consumed will be less 
than some fixed amount and the solenoid will give a fixed pull. 



Engineering handbooks can supply the following set of relationships: 



Pull 


IN 


(1) where N=total number of turns 


Current 


-j 


(2) where V=voltage applied 


Power 


~? 


(3) 


Resistance 


R = 4<< 2 /*L 


(4) where/ 7 = resistivity of wire 


Wire Length 


L = TT DN 


(5) 


Solenoid Diameter 


D = d+t 


(6) 


Coil Thickness 


f£ 


(7) where o< reciprocal of wire 

diameter (wire size) 


Number of Layers 


N 


(8) 



This is often the way a problem is presented to an engineer . If use of a 
computer is contemplated, the usual rules of analysis apply. That is, it 
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is best to do a little algebra before plugging in numbers . Too often the 
tendency is to choose wire size<< , start at equation 8 and work back, 
equation by equation, to arrive at power finally by equation 3. Applica- 
tion of algebra, however, gives: 

FV 

P = .-— 



4F/?je 



-<* 



!d ) 



(9) 



or power directly in terms of wire size. Savings in computer time can be 
just as important as savings in engineering analysis time. So it always 
pays to do a little analysis before muddying up the water with numbers . 

The important aspect of this problem is not that one can now compute the 
power output for a given size and thus choose an acceptable wire size, 
but rather that a computer program can be generalized such that, given 
any requirements for a solenoid, many different combinations may be 
tried quickly to find the suitable combination for the job at hand. Thus, 
to design a solenoid for a given pull F, it might be desirable to varyX , d, 
and V or, perhaps, to limit t. Using equations 1 to 8, a variety of opti- 
mizing programs can be written for solenoid design. The question which 
would determine the general form of each program is: What particular 
specifications are to be met? 



Heat Transfer Problem 



An insulating wall is made of three parallel layers of different insulating 
materials. The outside temperatures, t 1 and t , are known. The con- 
ductivity k. of each layer , a straight line function of the mean temperature 
t. , is known. The problem is to find the quantity of heat Q passing through 
a unit area of the wall per unit of time. 



ti 



*2 



*3 



t4 



x 2 



« - x 3 



Q is given by 



Q 



t i " h 



A-i An 



2 3 



(1) 
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However, the k. are dependent on the t. in the following manner: 



k-. = a-,t, + b. 



\ =3 2\ + h < 2 > 



k 3 = Vs + b 3 



where the a. and bj are given constants for the particular intervals . At 
the same time , the t- are given by 



vv(v 2 H (ti+t2) 

^t r (VS)=i(V'») 

2 2 

t n =t n -(^' t 4)=~( t 3 + M 



(3) 



but t 2 and t„ are not known. 

If the k. are known, t and t may be determined from the steady-state 
requirement that the quantity of heat passing through each layer per unit 
of time must be the same. 

Or 

k l (h ~ h) = k 2 <*2 " t 3> = k 3 $3 ~ W 

which gives 



(4) 







X 3 
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t 4 ) 
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x 3 
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+ — 


4- — -,,.. 
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k 3 
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= v 


x l 
k l 
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t 4 ) 


fi 


+ fi 


t f8 
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k„ 
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1 


2 


3 



But this produces a vicious circle; the k. are needed to determine t^ and 
to and vice versa. Thus, iteration is required. 

The procedure for this is as follows: 

1. Make a reasonable guess at t2 and t„. 

2. Use (3) to determine the t^. 

3. Use (2) to determine the kj. 
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4. Use (1) to determine Q. 

5. Since values are known for k^, recompute t 2 and t3 using equations (4), 

6. Repeat steps 2 through 5 until two successive values for Q agree. 

This problem converges to a solution nicely using the simplest iterative 
scheme, x = f(x,y,z). Experience has shown that the initial estimates of 
the internal temperatures are in no way critical except that t^ must be 
greater than t 2> which must be greater than t3, which must be greater 
than t4 . 

The data used in this sample computation is as follows: 



Layer 1 



Kaolin Insulating Brick 



b i 



.0623 



a. 
l 



.00010 



Layer 2 



Kaolin Insulating Firebrick . 0255 



.00005 



Layer 3 



Magnesite 



2.4395 



.00060 



1400 ( 



V 


= 200 u 


x l 


= 2.15" 


x 2 


= 2.0" 



6.0' 



The first set of results in Figure 6 shows the iterated solution to this 
problem giving the final value for the quantity of heat transmitted in a 
unittime as_Q = 25.38 Btu/min with the corresponding mean temperature 
t^, t 2 , and t3« An absurd case is also shown by the second set of results 
in Figure 6, with the input data changed in order to illustrate that con- 
vergence is not always possible. The data changed is 

H 
.02 

.01 

-.05 

Actually it can be seen that the reason for the divergence of the second 
se4^was^that4he^large^negative slope a3 c aused k3 to become negative (an 



impossible situation physically, since the coefficients of conductivity may 
not be negative). The value of -.05 for a 3 is not a reasonable value. It 
was chosen here to indicate that a mathematical procedure of this type is 
valid only within some range of selection of data. 
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PQ=0.0 

1 READ 10, Bl ,B2,B3,T1 ,T4,X1,X2,X3 

2 READ 10,A1,A2,A3 

3 T2=2.*(Tl-T4)/3.+T4 

4 T3= (Tl-T4)/3.+T4 

5 T1B=.5*(T1+T2) 

6 T2B=.5*(T2+T3) 

7 T3B=.5*(T3+T4) 

8 Z1=A1*T1B+B1 

9 Z2=A2*T2B+B2 
19 Z3=A3*T3B+B3 

11 Q=(T1-T40/(X1/Z1+X2/Z2+X3/Z3) 

12 PRINT 20,Q,T1B,T2B,T3B 

13 IF(ABSF(PQ-Q)-. 00005)2, 2, 14 

14 PQ = Q 

15 T3=T4+X3/Z3*Q 

16 T2=T1-X1/Z1*Q 

17 GO TO 5 

18 END 

Iteration Q ti to t 



3 



1 26.925558 1200.0000 800.00000 400.00000 

2 25.146838 1241.2234 671.36960 230.14618 

3 25.385972 1254.9914 684.25925 229.26788 

4 25.384732 1254.6856 684.23785 229.55224 

5 25.384376 1254.6691 684.21790 229.54884 

6 25.384383 1254.6698 684.21825 229.54845 



Iteration Q t-. t 2 to 



1 -384754.98 1200.0000 800.00000 400.00000 

2 188995.29 18589.197 83719.950 65930.755 

3 145.56898 853.61840 81.496675 27.878275 

4 38318.358 1390.8673 1208.5344 617.66710 

5 -823.77876 -77.502100-4718.9752 -3841.4731 

6 4076.2550 804.76095 192.05571 187.29476 

Figure 6. 
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Rocker Arm Cam Problem 



The accompanying figure shows a rocker arm and cam to be used in a 
fuel pump. 




Specified are: 

1. Axis of rotation of the off-center circular cam (X2, Y2). 

2. Eccentricity, E, and radius, R, of cam. 

3 . Axis of rotation of rocker arm (X.. , Y ) . 

4. Required angular rotation, oC , of rocker arm to give desired travel 
of pull rod. 

To be determined are: 



1. Necessary angle, B, between rocker arm and vertical. 

2. Necessary offset, L, for rocker arm. 

To obtai^a s^I^ 

expressions for B and L. For the "high" position of the rocker arm 

B = 90° - tf - arcsin ( L+R+E ) (1) 

A 

For the "low" position of the rocker arm 

L = A sin (90° -B- tf -</L )-(R-E) (2) 
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Note that because of the transcendental functions these equations may not 
be solved explicitly for either L or B. 

The direct iterative procedure for solution is: 

1. Make a reasonable guess at L. 

2. Determine B from equation (1). 

3. Then use this value of B in (2) to obtain a value for L. 

4. Repeat steps 2 and 3 until two successive values for L agree. 

This model requires iteration for the basic variables B and L. In actual 
practice, many other parameters are to be specified. These can be com- 
puted directly once B and L have been determined. 

The solution here is interesting because the straightforward iteration 
procedure diverges. If equations (1) and (2) are stated as: 

B=f!(L) 

L = f 2 (B) 

then a successful iterative equation to replace equation (2) in step 3 above 
is: 



L t+1 = L i - f 2( ] ^L+l) - L i (3) 



where 



B i + i =f i< L i> 

Using this iterative procedure, convergence is relatively slow. An 
attempt to improve the convergence using the Newton-Raphson technique 
gives a divergent iterative scheme. However, a modification of the pro- 
cedure shown does improve convergence. 

The improvement was made simply by changing the equation to read: 

L i+1 = L i " ( f 2( B i + l> " L i> < 4 > 

This is treading on obviously dangerous ground, since arbitrary changes 
in the iterative equation are not easily justified. For this problem, how- 
ever, it does work. This is mentioned here to show that experimentation 
with techniques is necessary in some instances when other approaches 
fail. But keep in mind that it is important to have some way of judging 
the correctness of the results. 

Consideration of the following graph may assist in understanding the 
iterative process. Let 

x n+l = * ( x n) 
be the direct iteration equation, and represent all other schemes as 
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x n+1 = W(f (x n ), xj 

where the W represents the weighting function used in making the next 
estimate. 

In either case the evaluation of the function giving the next estimate might 
be plotted as in the graph. 

For a solution, the next estimate of x and the input value of x must be 
equal ( f (x) and x as shown in the graph). Assume an error, € , is made 
in the first estimate, x . The graph indicates the next estimate to be in 
error by an amount C £ where fc 2 > € 1 (^ e technique is diverging and 
€ 3 is even greater than € 2) • 

Given the initial estimate x n and the graph as shown, an obvious proce- 
dure for preventing the divergence is to make use of the difference 
( 6 9 - 6 ..) in writing a new iteration equation — that is: 

x n + l =x n-< 6 2- 6 l> 
and since ( 6 2 _ € ■]) = f ( x n ) - x n 
then 



x .1 = x " (f ( x ) ~ x ) 
n+1 n v x n' n' 




61 (62-ei) (€ 3 -€ 2 ) 



This is recognized as equation (4) above for the cam problem. If 
(6 2 ~ € ,)> €.. this procedure overshoots the correct result and 
divergence may again result. Equation (3) in the cam problem used a 
weighting factor of 1/2 to prevent the overshoot. Choice of the "best" 
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weighting factor would require complete knowledge of a similar graph — 
which of course is not available without a solution. 

The intent here is to indicate the reasoning which may go into the selec- 
tion of iteration equations. 

The direct and Newton-Raphson iterative techniques are not the only pro- 
cedures available. For example, wide use has been made of a procedure 
credited to Mr. J. H. Wegstein of the National Bureau of Standards. The 
Wegstein method does not require evaluation of the first derivative and 
therefore is a powerful tool in cases where the first derivative is im- 
possible or difficult to evaluate. 

The rocker arm cam problem illustrates the use of experimentation to 
find a successful iterative scheme . No mention has been made of the use 
of an iterative technique to handle simultaneous algebraic equations . 

The problems discussed in this manual have been treated in the form 

x = f (x, y, z) 
Physical problems may frequently give rise to forms similar to: 

x = f x (x, y, z) 

y = f 2 (x, y, z) 

As a matter of fact the rocker arm cam problem can be treated as a prob- 
lem giving rise to simultaneous algebraic equations. 

The Newton-Raphson and Wegstein methods may both be extended to han- 
dle simultaneous equations. 



Spring Problem 



The engineering design of springs makes an excellent computer applica- 
tion, especially from the economic viewpoint. It takes a great many 
springs to keep this mechanical world operating, and a less expensive or 
longer-lasting spring design can mean large returns dollar-wise. This 
example is concerned with the design of a particular type of spring, but 
the general method of computer solution would apply to almost any spring 
problem . 
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i.DDi. 



DD 



Cross Section of a 
Double Helical 
Spring 



TD DT 



DD 



A double helical spring is to be designed to support a given torque winding, 
such that the mean diameter and length of the spring will be within set 
limits — the shorter the length the better. The variables specified (sub- 
script 1 refers to inner spring, subscript 2 to outer spring) are: 

M-l , M 2 = Winding torque 

S 1 , S ? = Maximum selected stress, or ultimate stress, of material 

T^ , T 2 = Maximum torsion angle 

K.. , K 2 = Required spacing between turns 

C - Required spacing between coils 

E = Young's Modulus 

The following are to be computed: 

h-^ , h 2 = Dimensions shown in the diagram 

b-i , b 2 = Dimensions shown in the diagram 

D l ' D 2 = Diameter of springs 



■^1 ' L 2 = Length of springs 

Njl , N 2 = Number of turns in the springs 

Known relationships are: 
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(1) M = E bh T Gives winding torque for a given torsion 

6-6DN ang ie. 

Sbh2 

(2) M = — - — Limits winding torque in terms of ultimate 

stress of the material. 

(3) Lg = (N 2 + l)(b 2 + K„) Outer length by geometry. 

(4) D;l = D 2 -h^-h2-2C Inner dimension by geometry. 

Note that for the outer spring five variables (h 2 , bg, D 2 , Lg, Ng) must 
be determined but there are only three equations which apply. One 
approach to obtaining a solution is: 

(a) Set Dg (because limit exists) 

Ng (because this must be either an integer or 

an integer plus . 5 — successive guesses 
will be simpler) 

(b) Equate (1) and (2) and solve for h. 

= 1. 1D 2 S 2 N 2 
2 ET 2 

(c) Then from (2) obtain 

6M 

b 2 = S < 6 > 

Sghg 2 

(d) and Lg = (Ng + l)(b 2 + Kg) (7) 

This provides the parameters for the outer spring. At this point it is 
necessary to make an additional guess to handle the inner spring. 



(e) 


Set N-l 


(f) 


Then solve 




D 1 = Dg-h 1 -hg-2C 


and 


[ h^i-iDlSlNl 




ET X 



simultaneously for h^ obtaining 

h = 1. 1 S 1 N 1 (D 2 -h 2 -2C> 
1 ETj^ + i.iSjNj (8) 



(g) From (2): 



6M-, 
b = S (9) 

1 s lhl 2 
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(h) and finally: 

L x = (N 2 + l)(b 2 + K 2 ) 

Now testing is performed to see whether the springs are suitable; L 2 
must be equal to L-. 

According to the design criterion — that L 2 approximately equal L-. — a 
decision can be made to accept or reject this design. In case of rejection 
the selection for D 2 , N 2 and N]_ can be modified and the indicated cal- 
culations repeated. 

A logic diagram of the procedure is shown in Figure 7. 



Set 
D 2 , N 2 , Nj 



o 



For Outer 
Spring Calculate 

h 2' b 2' L 2 



For Inner 
Spring Calculate 



Select 

New 
D 2 , N 2 , N x 




Satisfactory 
Design 



Figure 7. 



The question of how to modify the selection of D 2 , N 2 and N-^ must be con- 
sidered. The mathematical nature of the problem (nonlinear equations 
with an infinite number of solutions) does not permit use of the iterative 
techniques mentioned thus far. Two possible procedures for varying the 



selection may be called scattering and search. Both of these procedures 
assume the existence of a maximum or minimum for the function G which 
is used to denote "goodness" or "acceptability." In the spring problem 
the function is 



G- 



L 2" 



L x I = G (D 2 , N 2 , N x ) 
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Considering G as a function of only two variables 

G = G (D 2 , N 2 ) 

allows visualization of the graph of G shown here. A corresponding three- 
dimensional picture for 

G (D 2 , N 2 , N x ) 

could also be visualized. 

The continuous curves represent contours on which G is everywhere equal. 
There may exist one or more minimums as shown. 

Scattering 

A set of selected values for D 2 and N 2 (in three dimensions, D 2 , N 2 , and 
N-i ) which create a grid covering the possible range in both D 2 and N 2 is 



T 



Allowable 

Range for 

D_ 



O 1st trial points 
A 2nd trial points 




Allowable Range For 

No 



Figure 8. Scattering 

used to compute corresponding values for G. The grid is shown in Figure 
8 as a set ofO's. A smaller area is selected around the minimum G and 
the process is repeated, as indicated by the A's in the figure. 

This corresponds to the familiar trial-and-error approach except that 
many tries may be run at one time on a computer. Each solution may be 
made available allowing consideration of secondary design criteria. In 
the spring problem, a minimum | L 2 - L | is not sufficient for a good 
design; the ratios j^ > _h£ anc ^ ^1 are a l so important. It is quite probable 
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that the finally selected spring would have a small but not the smallest 



L 2" 



and other desirable features, 



Search 

One starting point in n dimensions is selected (shown in Figure 9 for the 
two-dimensional case). 




Figure 9. 
Search 



Start 



For each variable in turn: 

(a) The variable is changed by some small amount and G is recomputed. 

(b) If G decreases, (a) is repeated until such time as G begins to increase 
and then the variable is set to correspond to the smallest G. 

(c) If G increases for (a) , an equal change is made to move in the opposite 
direction in a manner similar to (b) . 

Steps (a), (b) and (c) are repeated for each variable in turn, until such 
time as new changes within the practical bounds fail to decrease G sig- 
nificantly. 

This method moves more quickly to the solution and represents a more 
fully automated procedure. However, it demands that the behavior of the 
function (in this case G) be well understood. Discontinuities or lesser 
"minimums" can destroy the effectiveness of this method. 



The scattf 



">ften be 



sJbehavioi 



criterion function such as G (D£, N£, N-^) above, to the extent that the 
search technique may be used in succeeding computations . 

It is possible to add other selection criteria in the search technique. For 
example, the material cost in designing the above spring might be used 
such that the selected spring has the lowest material cost for those designs 
with a value of G within a specified range . 
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Pitman Arm Stress Analysis 



A problem often encountered in mechanical engineering is that of design- 
ing a member to withstand given tensile and shear stresses. This partic- 
ular example concerns a pitman arm (see Figure 10) as a part of an auto 
steering mechanism. The drawing shows a pitman arm with the applied 




SECTION AA 



Figure 10. Pitman Arm 

force F which gives rise to the various stresses in the arm. 

The usual approach to this problem is for the engineer to specify dimen- 
sions which seem appropriate and then to check the design with calcula- 
tions. This is done by choosing cross sections of the design at intervals 
along the arm — computing the stresses at each cross section and check- 
ing that these stresses are less than the failure stress for the material 
used. 

For Section AA shown, the computations are: 

(a) For bending stress 
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3 b ='£k (1) 

Z x 



For torsional shear stress 

z 

p 



where Z x and Z are called section moduli and depend upon the moments 
of inertia for the member. For the elliptical case the moduli are 



z x =f 2 bh 2 (3) 

Zp - f 6 bh 2 (4) 

(b) The total normal tensile stress and the total shear stress can be com- 
puted in terms of the bending stress and the torsional shear stress. 

For normal tensile stress 

L\ 2 2 1/21 

S n = 2 LSb +<S b + 4S t ) J (5) 

For shear stress 

12 2 1/2 

S s = ~2< S b +4S t> < 6 > 

Once these total or maximum stresses have been calculated, they are 
compared with acceptable values for the material to check whether a 
feasible design has been arrived at. 

An approach which may be of more use to the engineer — and which is 
facilitated by the use of computers — is to start with the maximum allow- 
able stresses and from this to obtain the dimensions required for the 
member to support these stresses. 

For example, assume that the pitman arm shown is to be elliptical in 
cross section, with h=2b, and that values have been assigned to S n and 
S g . The problem now consists of determining values for b and h in terms 
of these values. 

Equations (5) and (6) can be manipulated to give S b and S^. in terms of S n 
and S s : 

Sb = 2(S n - S s ) (7) 



T 2 2l 1 /2 

S t = L S s - (S n -S S ) J (8) 

Equations (1) to (4) may be used finally to restrict b and h in terms of 
S b and S-: 
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and 



" 9 

— ui — 






and 



b d -_8 Fl 2 = 
b 3 -_4 Fix = 



or, since h = 2b, 

(9) 
(10) 



This gives two cubics which restrict b. The proper design must have a 
value of b which is at least as large as the larger of the two roots obtained 
from equations (9) and (10). 



Now the proper dimensions can be computed for the arm as we move along 
the arm (that is, change 1-^ and 1 ? ). Furthermore, if the curve used for 
the arm can be expressed mathematically, then 1^ and 1q can be auto- 
matically computed as the arm is traversed toward or away from the axis 
of rotation. If this is impossible, then 1-. and 1 2 must be measured for 
various cross sections of the arm and the corresponding dimensions com- 
puted for the cross sections. 

The point of this example is that very often a computer approach should 
start with a reanalysis of what is to be done, considering the possibilities 
of improving the usefulness of the results. 



Automobile Suspension 



The following three-dimensional problem occurs in the engineering of auto 
suspension systems; arms of length d^, dg and d 3 are fixed at points in 
space A,B and C respectively. 



(x,y,z) 
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The point of juncture of the three rigid arms is found to be the intersection 
of three spheres with centers at A,B and C and radii d-j_, d 2 and d3 respec- 
tively. 

Exprp?sing this in equation form gives the following three equations: 

2 2 2 2 

(x - x x ) + (y - y x ) + (z - Zl ) = d]_ (1) 

(x - x 2 ) 2 + (y - y 2 ) 2 + (z - z 2 ) 2 = d 2 2 (2) 

(x - x 3 ) 2 + (y - y 3 ) 2 + (z - z 3 ) 2 = dg 2 (3) 

This set of equations appears difficult to solve for x, y and z, and iteration 
has been used in many instances to find a solution. Iteration in this prob- 
lem is time-consuming and unreliable (convergence is not always possible). 

If the following algebraic steps are performed, an explicit relationship for 
x, y and z is found: 

(a) Perform indicated squaring in equations (1), (2) and (3). 

(b) Subtract equation (2) from (1) and similarly equation (3) from (2) . 
This will eliminate terms in x 2 , y 2 , and z 2 giving two linear equations 
in x, y and z. 

(c) Solve these two linear equations algebraically to obtain x and y as 
linear functions of z. 

(d) Substitute these functions for x and y into equation (1) and solve the 
resulting quadratic in z. 

This will then give equations of the form: 

9 .1 

z = -b + (b z - 4ac) 2 

2a 

x = e + f z 

y = g + hz 

where a,b,c,e,f,g ana n are intermediate substitutions for terms made 
up of the known quantities xj, yj, z^ and dj. 

This gives two complete sets of answers: the correct physical answer and 
its mirror image as illustrated by the dotted lines . The selection of the 
proper set of results can be accomplished by programming the logical 
considerations of the physical systems 

A computer is convenient for the evaluation of the lengthy equations in- 
volved. If the problem is to accomplish this evaluation for many points 
in the motion of one of the points (say point C, as indicated by the arrow), 
the use of a computer becomes mandatory. The orientation of a suspen- 
sion system must be known for the many possibilities in the auto's motion. 
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This problem illustrates an important class of mechanical engineering 
problem — that of determining the orientation of physical bodies by speci- 
fying the coordinates of points in the body. 



Spring Stress Analysis Problem 



The following illustrates the use of the same algebraic approach in another 
problem which also involves the use of the summation of forces so familiar 
in stress analysis. Find the necessary tension, S, in the spring for the 
support of W at various values of the angle -9-. 




The analysis might proceed as follows: 
To find point A 



x l ~ ^1 cos "®* 



y. = d„ sin -0- 



To find point C 



(x-x-l) + (y-y-L) = d 3 
(x-x 2 ) 2 + (y-y 2 ) 2 =d 2 2 



(1) 
(2) 

(3) 
(4) 



Following the procedure as outlined in the previous problem, expand (3) 
and (4) to give: 



x - 2XXJL +x 1 + y - 2yy;L + y 1 = d 3 
x 2 - 2xx 2 + x 2 2 + y 2 - 2yy 2 + y 2 2 = d g 2 
Subtract (6) from (5) to give: 



(5) 
(6) 
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2x (x 2 - Xl ) + ( Xl 2 - x 2 2 ) + 2y (y 2 -yi) + (yi 2 -y 2 2 ) = ( d 3 2 " d 2 2 ) 
or x = Ey+F (7) 



Where 



E 



frry2> 

< x 2" x l ) 



and F = [(d 3 2 -d 2 2 ) - (y^-y^) - (x^-Xg 2 )] /2 (x 2 - Xl ) 

Substituting from (7) into (5) gives 

E 2 y 2 + 2EFy + F 2 -2x 1 Ey-2x 1 F + x 1 2 +y 2 -2yy 1 -y 1 2 = dg 2 or 
(1+E 2 )y 2 + (2EF-2x 1 E-2y 1 )y + (F^fcKjF+x^-y^-dg 2 ) = 
or finally y = - Q + Vq 2 -4PR 



2P 



where 



P = (1+E 2 ) 

Q = (2EF - 2x 1 E-2y 1 ) 

R = (F 2 -2x 1 F+x 1 2 -y 1 2 -d 3 2 ) 
and x = Ey+F 
To find angle (J) 

x„ = (d 1 + d.) cos -0- 

y 3 = (d x + d 4 ) sin -0- 

and (J) = tan" (y-y 3 ) 

The determination of (J) allows the use of the force diagram and the equi- 
librium principle, which states that the sum of forces acting in both the x 
and y directions must be equal to zero, to determine the necessary tension, 
S, in the spring. 




'W 
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or 



£ F x = T cos-9- - S sin ( (J) - 90°) = 

£ F = T sin-0- + S cos ( - 90°) - W = 

T = S sin ( (J) - 90°) 
cos-0- 

S = W/ [cos ( (J) -90°) + tan -0- sin ( (J) -90°)] 

The usefulness in the design of a program to evaluate the tension in the 
spring for any configuration, motion of-0- , or load W is obvious. In this 
case the evaluation depends upon locating points in a plane. Traditionally 
this type of determination has been made on a drafting board . As in this 
case, where motion or many possible designs are to be evaluated, cal- 
culations on a computer can greatly reduce the time invested. The improve- 
ment in accuracy can be even more important. For example, a one-degree 
error in the steering geometry of an auto can completely change its steering 
characteristics. 



Vehicle Simulation Model 



It would be unfair to close this brief discussion of a few engineering 
mathematical models without mention of at least one of the larger models 
which have been constructed. 

Several auto manufacturers have succeeded in representing the operation 
of cars and trucks by a computer model . By operating the model in the 
computer, such items of performance as acceleration, fuel economy, and 
riding characteristics can be determined for a proposed model before the 
model is built . 

The procedure for accomplishing this introduces a very interesting con- 
cept. Normally the solution to a set of equations is obtained by solving 
for one point in time — a static situation. The performance of a car, 
however, is a dynamic situation. Beyond this, equations to describe all 
the myriad operations in moving a car simply do not exist. Rather, there 
exist combinations of theoretical equations , empirical equations and tab- 
ular data to describe the functioning of various parts of a car. The heart 
of a car simulation model is the logic which ties these pieces of theory and 
data together under one unifying variable . Different approaches choose 
different unifying variables; for illustration, however, suppose the variable 
chosen is time . 

Typically, the motion of the car might be studied time-wise for a set level 
of fuel input. For the first unit of time the information representing the 
relationships between parts of the car is interrogated for the motion to be 
expected, working from the engine back through the transmission and 
finally to the wheels . Once the motion for the first unit of time is deter- 
mined, the time is incremented and the operation repeated. 
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Another approach might be essentially the reverse of this: Starting with 
the desired velocity of the vehicle, work back through the wheels, dif- 
ferential and transmission to determine what the engine speed must be to 
give this velocity. However, time is again the unifying factor in that con- 
ditions for the previous time interval help determine the required motor 
speed for a particular velocity. 

A look at the basic relationships in this example may be helpful . 

Assume we want to know the linear acceleration of the vehicle at full 
throttle as a function of time . 




Starting with a given engine speed N^ and linear velocity V, the first set 
of relationships is concerned with the engine, A. 

(a) The first step is to use a table to determine the output torque of the 
engine . 

T i = f i < N i> 

The torque delivered to the transmission may likewise be a tabular 
function of N^: 

T 3 = f 2 (Nj) 

as is the torque consumed in friction of the engine parts: 

T 2 = f 3 ( N l) 

At this point the acceleration of the engine for this time interval may 
be determined as 

T r T 2- T 3 
* m = — 

where I = total engine inertia 



(b) Drive shaft speed, N„> is computed by 



N 2 =.V_ G 



where G = axle ratio 
and L = rolling radius . 
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The transmission ratio, TR, is a function of N 2 /INL 

TR = f (N 2 /N 1 ) 

This relationship is best represented by a table of TR vs N«/N 1 which 
has been determined experimentally for the particular transmission 
to be used. 

(c) Moving on to the drive shaft, the torque delivered to the drive shaft 
may be computed as 

T 4 = T 3 • TR 

(d) Finally, arriving at the rear wheels, the force exerted on the wheels 
due to the delivered torque may be given as 



G_ 

L 



F =— <T 4 - T 5 ) 



where T 5 = torque loss due to friction and acceleration of drive line 
parts . 

Then the acceleration of the car, a, is given by Newton's second law 
of motion 

a = F - R 

M 

where M = mass of the car and 
R = rolling resistance. 

The rolling resistance is generally given as an empirical function 
such as 

R = aM + b AV 2 

where A = frontal area of vehicle, and a and b are experimentally 
determined constants (see section on empirical relationships) . 

At this point the time is incremented and a new engine speed N- is com- 
puted as 

^m = *m + A * m+1 



Q*i> ,r( N iL + ^ At 

1 m + 1 J- m \ m 



m 



The above procedure A,B, C and D may be repeated for as many time 
increments as desired. 

This logical procedure may be represented by the flow chart in Figure 11. 
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The flow chart in Figure 12 is included to illustrate how essentially the 
same mathematical model may be used in reverse to compute fuel economy. 
In this case the desired motion as a function of time is specified. Cal- 
culations determine the engine speed required to give the desired motion 
for each time interval. Finally the fuel consumed at a given engine speed 
is known and therefore the sum of the fuel to be consumed for all time 
intervals may be calculated. 

Many simplifications of vehicle simulation have been made in this descrip- 
tion. 

The major objective in presenting a description of vehicle performance in 
this chapte r is to illustra te the use of the three types of information 
(theoretical equations , empirical equations , and tabular test data) under 
a logic format to obtain a computer simulation of the physical object. 

The next chapter describes the use of the computer as a tool in finding 
empirical relationships . Actual procedures for handling of tabular data 
in a computer program are also discussed. 
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Exercises 



1. Perform the algebra required to obtain equation 9 in the solenoid de- 
sign problem , using equations 1 to 8 . 

2. Perform the algebra indicated in the suspension analysis problem and 
either show the complete equations for x, y and z or show the equa- 
tions which define the intermediate variables a, b, c, r 5 s, t and u. 

*3. Numerical techniques are often used to solve simultaneous differential 
equations with nonlinear coefficients . Prepare a logic flow chart to be 
used to program the solution of: 

d2 y ,dx 

J +oC — +xy = A 



dt2 dt 

dx+dy +<3 < =B 
dt dt 

where A and B are known constants and oC - t + sin << t , given t Q , x , 
y and fdy\ . 
\dtj 

Hint: Use Euler's approximations as described in problems I and J 
of Section IV for the differential equations and Newton's procedure 
for evaluating &£ at each interval in t. 
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SECTION VI: SELECTED MATHEMATICAL TOPICS 



The problems discussed in this section have been chosen to answer some 
of the legitimate questions that many practicing engineers might well ask 
about terms often heard in computer groups . Such things as matrix inver- 
sion, relaxation techniques and eigenvalues are out of the realm of expe- 
rience of many engineers. It is felt that some simple examples illustrat- 
ing these terms may lead to eventual recognition of real problems which 
can be tackled using the techniques described. 



Matrix Inversion 



A major application of computers is in handling the solution of large sets 
of simultaneous equations which may occur in such engineering areas as 
stress analysis, statistical least squares and circuit analysis. 



1 f ^ *3 ] \ 12 $ R 2 



For the simple circuit shown, applying Kirchhoff's law gives a set of linear 
equations in i and i : 

J- £t 

(Rj + R 3 ) ^ - R 3 i 2 = Ej 
-Rgi! + (R 2 +R 3 ) h = - E 2 

This set of two equations in two unknowns does not require sophisticated 
handling. However, consider a slightly more complex circuit in which 
the values of the resistances are known and the currents are to be deter- 
mined. Again, Kirchhoff's law may be used to establish a set of linear 
equations: 



Ii +I 9 +I« 



W^IO 



~h + h ~ h 



-Iq + U ~ I( 



hl 1 ^^ 



1 amp 



-R 7 I 7 + R 8 I 8 - R^io - ° 



-R 5 I 5 + R 8 I 8 - R 9 I 9 - 




2 amps 



R 2 I 2 " R 3 I 3 " R 5 J 5 - ° 
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-R-^ + R 2 I 2 - R 4 I 4 



-R4I4 ~ RqIq + R7I7 



Since the solution of this set would require a considerable expense of time 
and effort, the problem is well suited for a computer. 

The general procedure for solving a set of simultaneous equations on a 
computer is to make use of matrix inversion techniques. A standard 
elimination technique for matrix inversion is performed in the following 
manner: 



Given the coefficient matrix 

a ll a 12 a 13 • • • a ln 
a 21 a 22 a 23 • * * a 2n 



a nl a n2 a n3 ' ' ' a nn ° 



□ □□...□ 



add a column — a unit vector — which contains a 1 in the first row and 
zeros elsewhere. At the same time, add a row — called the pivot row - 
denoted by CZI ' s . Then perform the following computations to arrive at 
a new array: 

1. For the pivot row elements: 



a . = a l, j+1 
P,J — — 



j = 1,2 n 



Xi 



2. For all other elements, compute a new value: 



= a 



i»] 



i,j+l 



- (a. J (a .) i= l,2,...,n-l 



i>l P,j 



j= 1,2, ...,n 



3. As a result of step 2, all the new elements of row 1 are zero. This 
row is dropped, and the remaining n rows renumbered 1 through n. 
Thus, for the last row: 



a' = a 
n,j p,j 



j = 1,2, ...,n 



4. Add a new unit vector and pivot row and repeat steps 1, 2 and 3 a total 
of n times. The resulting array is the inverse of the original matrix. 

Given a set of simultaneous equations: 
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a ll x l + a 12 x 2 + . . . + a ln x n = b x 



a - Xl +a n2 x 2 + ' ' ' + a nn x n 



"nl~l 



then the solution can be obtained directly by starting with an n+1 by n array 
in which the original matrix is augmented by the b vector. If values of x 
are required for more than one set of b values , the additional b vectors 
can be incorporated in the original array as shown: 

a ll a l2 • • • a ln b ll b 12 
a 21 a 22 • • * a 2n b 21 b 22 



a n l a n 2 • • • ann b n i b n 2 

Subjecting this m by n matrix to the four-step procedure above a total of 
m times gives the array: 



x ll x 12 a ii a i 2 



x x a 1 a' ... a' 
21 22 21 22 2n 




\ . 

\ x „ x n a' a' . . . a' / 
\nl n2 nl n2 nny 

where the a! . are the elements of the inverse of the original coefficient 

matrix, and the x^ j are the solutions for each of the two b vectors of the 
matrix equation 



AX =B 



Eigenvalue Problem 



This problem is discussed to show, with a simple example, an application 
of computers in handling sets of differential equations . It illustrates the 
meaning of eigenvalues for a set of differential equations which describe 
the motion of a mechanical system. 



4 



x l 



'/ — dSULWb—C^) — <mHr^\-dMMs^^)-^mJlM^- 



x 2 



Tn. 1 = m.r, : =in.' 



/ *"l-*"2-" 1 3 
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Given the spring and mass system shown, assume that the normal modes 
of oscillation of the system are to be determined. 

The differential equations of motion can be written: 

m^x-L = - kx-L + k(x 2 - x^) = k(x 2 - 2x^) 

m 2 x 2 = - k(x 2 - x x ) - k(x 2 - x 3 ) = - k(2x 2 - x x - x 3 ) 

m^Xg = - kx 3 + k(x 2 - x 3 ) = k(x 2 - 2x 3 ) (1) 

A procedure for solving these differential equations is to assume solutions 
of the form 



Xl = Ae lwt , x 2 = Be 107t , x 3 = Ce iwt 



(2) 



Substituting these into equations (1) , performing the indicated differentia- 
tions and rearranging terms gives: 



(2^-w 2 ) A 



k R 
m B 



~~ A + (2|-a J 2 )B- ^C 



= 

= 



(3) 



m" B + (^ - w 2 )C = 



v m 



which can be written in matrix form as 




or more simply, 

(D- M)X=0 



(5) 



(4) 



For A, B, C to satisfy equations (3), the determinant of the coefficient 
matrix must vanish. This is equivalent to the statement 



(D- M)=0 



(6) 



Evaluating the determinant for equations (3) and equating it to zero gives 
a polynomial — called the characteristic equation — in oj 2 : 

2 
# - oj 2 ) 3 - 2^ (2^ - UJ 2 ) = (7) 



x m 
with roots 



m v m 



2 2k 

1 m 

2 2k 



a) 



= — + V 2 



2 m 



m 



(8) 
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The values of w which satisfy these equations are called eigenvalues. In 
general, values of X which satisfy (5) are eigenvalues. The vector X is 
called the eigenvector. 

It can be seen that for larger sets of differential equations the solution of 
the characteristic equation requires some means for computing the roots 
of a polynomial. 

For this problem the eigenvalues give the natural modes of vibration for 
the mass and spring system. This calculation of eigenvalues for differen- 
tial equations is termed frequency analysis. In an earlier section attention 
was given to the solution of differential equations . In motion problems 
this amounts to determining the displacements as a function of time and 
is termed an amplitude analysis. Both frequency analysis and amplitude 
analysis are important computer applications. 

To illustrate a computer solution to a frequency analysis , consider the 
spring and mass system shown in the diagram. 




The differential equations which describe this system are as follows: 
- + a ll x l + a 12 x 2 + a 13^1 + *L4?2 = ° 



d x- 



dt 2 



A 2 

d xo 

- — +a x+a x+a -0- +a -9- = 
,.2 21 1 22 2 23 1 24 2 

dt 



2 
d "©--. -9- -9- 

~ + a 31 x l + a 32 x 2 + a 33 1 + a 34 2 
dt 2 



d 2 e . 

dt 2 



2 + a 41 X l + a 42 X 2 + a 43 l + a 44^2 = ° 



where the a. . are dependent on the spring constants and the masses, 
Assume solutions of the form: 



x^ = x, ^ COSOJt -9- = -9- cos cut 

1 10 1 10 
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x 2 = x 2Q cos cut 



*2 = ^20 COS ^ 



where x-^q, x 2 q, -0-iq and -0- 2 o are ^ ne initial displacements. The appro- 
priate differentiations and substitutions give a homogeneous set of linear 
equations of the form: 



a ll a 12 a 13 a 14 



a 21 a 22 a 23 a 24 



a 31 a 32 a 33 a 34 



a 41 a 42 a 43 a 44 





x 10 




X 20 




^lO 




" 9 '20 



CO 



k 10 



k 20 



10 



r 20 



The iterative procedure for the solution of these equations involves the 
following steps: 

1. Initially, let x 10 = x 2Q = -©- 10 = " e -20 = 1 -° and evaluate the left-hand 
side for new values of ui^x-^q, a> 2 X20, w^-B-.q and oj 2-0- 

That is: 

2 
cjx =ax +ax +a-9- + a -e- 
10 11 10 12 20 13 10 14 20 

and so on. 

2. "Normalize" for new guesses at x , x , -e- , and -0- 9O by setting 



x 10 1 



tu -e- 



10 



,2 
Cd X 



10 



CO X 



10 



X 20 



20 



2 

CO X 



10 



2 „ 
co -0- 



20 



k 20 2 

CO X. 



10 



3. Repeat steps 1 and 2 until such time as successive values of x-^q, x 2 „, 
-0- , and -e- are very close. At this time, convergence has occurred 



and o>2 can be computed. 



2 2 

4. Upon convergence oj can be evaluated from the given value of co x-,q\ 



say 



(J x = K 
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The following is the iterative solution for oj 2 of a sample set of coefficients: 
Particular equation set 



+929 


-332 


+3,590 


-748 


+748 


-8,090 


+4.17 


-4.17 


+35,400 


-14.31 


+14.31 


-40,300 



Iterative solution 



-3590 




x 1 




x 1 


+8090 
-11,760 




x 2 


= u> 2 


x 2 


+40,300 




- 2 . 




*2_ 



Iteration 


x l 


1 


1.0 


2 


597.0 


3 


143085.0 


4 


76883.0 


5 


63414.0 


6 


60681.0 


7 


60002.0 


8 


59825.0 


9 


59777.0 



1.0 
0.0 
-321095.0 
-171908.0 
-141558.0 
-135398.0 
-133869.0 
-133468.0 
-133362.0 



1.0 
23640.0 
1401773.0 
477975.0 
349238.0 
323511.0 
317130.0 
315460.0 
315017.0 



1.0 
0.0 
-1595813.0 
-844314. 
-693154.0 
-662485.0 
-654870.0 
-652876.0 
-652347.0 



oT = 59777.0 



Knowing that x in is to be set to 1, then 



u) 2 = K 



To clarify what has been done here, remember that the set of equations 
has no constant term — it is a homogeneous set of equations. Essentially 
this means that there are an infinite number of solutions which satisfy the 
equations. This is certainly reasonable when the physical system under 
consideration is examined. In a vibration problem of this kind the initial 



displacements x-^q . 
ferent values . 



L 20 = 



20' 



and B"2o must be expected to take on dif- 



In the above procedure a value of x-^ was selected which then fixed the 
values of the other variables x 2 q, "^io> ^20 an< ^ allowed a; to be deter- 
mined. 

A very similar iterative scheme (without the normalization step) can be 
used to solve sets of nonhomogeneous linear equations . 



Partial Differential Equations 



Many engineering problems involve the handling of partial differential 
equations. Three classes have been distinguished as: 

1. Elliptical equations (describing potential fields) 

2 
V (J) = g (x, y, z) 
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2. Parabolic equations (describing heat flow and diffusion) 

V 2 (D =kil 
it 



3 . Hyperbolic equations (describing wave action) 



V 2 <D 



at 2 



2 . 
if 

V~0 



where y is the Laplacian operator in rectilinear coordinates: 
2 - * 2 d> + ^20 ^20 



Jx^ 



^ 



*z' 



A basic approach to handling partial differential equations when describing 
a particular material or space is to create a grid of points covering the 
space as shown for a two-dimensional space . 



Then, at any point the 
first derivative with respect 
to x can be approximated in 
two ways: 



*<h ^s-f 



*x K 



Ax 



or 



(1) 



\§ 



\ *x / 



4>o" 4>: 
Ax 



2 

AX 




44 



The second derivative can be approximated as 



*x 2 



^ " <Ti> 2 ♦s + ^i-^ 



Ax 



(Ax)' 



(2) 



Derivatives in the y direction can be obtained in the same way. Using 
this procedure, any partial differential equation may be reduced to dif- 
ference equations amenable to computer handling. 



Potential Problem 



Apply this to the two-dimensional problem of finding the potential dis- 
tribution in a square whose sides are maintained at voltages 

< V 1>B' < V 2>B' ^B' ^ d < V 4>B 

If there is no charge within the square, the potential distribution is de- 
fined by the Laplace equation: 



* X 2 * }y2 



(1) 
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< V 3>B 



(V,) 



l'B 



(V ) 
V 2 ; B 



< V 4> B 



A 



Setting up a square grid system to cover the square , for the general 
point A, the approximations for the partial derivatives become: 






Va-V, 



OX) =v -v 3 

^x 2 



UL- V 4 + V 3 - 2V o 



}x 2 



Similarly, for the y dimension: 



u 2 



V + V - 2V^ 
J 1 o 

h 2 



Then (1) becomes: 



V l + V 2 + V 3 + V 4 " 4 V o = ° 



(2) 



This is the basic relaxation equation. It is applied in the following way: 



1. A first guess at the potential of each point on the grid is made on the 
basis of the known boundary conditions . 

2. Moving systematically through the points of the grid, the quantity 
called the residual is computed for each point and stored. The residual 
is given by: 
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Figure 13. 
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R o ^ V l + V 2 + V 3 + V 4 - 4V o 

Initially, equation (2) will not be satisfied since the potentials are 
guesses. 

3 . Again moving systematically and considering each point not on the 
boundary, the potential is adjusted to make the residual for the point 
equal to zero by applying the following equation: 

newV = V 0+ R /4 

4. This affects the residuals of the surrounding points, so they are ad- 
justed by: 

new Rj = Rj_ + R Q /4 

5. Steps 3 and 4 are repeated until no residual is found whose absolute 
value is greater than some predetermined limit of accuracy. At this 
time the relaxation equation is satisfied and the potential distribution 
is known. 

It is possible to quickly write a FORTRAN program to do the necessary 
computation. For this problem let 

M = Number of points in grid on x axis (200 max.) 

N = Number of points in grid on y axis (200 max. ) 

V (I, J) = Potential at points on grid (initial guesses) plus boundary 
values 

R (I, J) = Associated residual 

DEL = Limit of accuracy desired. 

The resulting FORTRAN program is shown in Figure 13 . 

With this example in mind, it is of value to consider how a more complex 
problem of this type might be handled. 



Temperature Distribution Problem 



For a cylinder made up of layers of material of different conductivities, 
find the steady-state temperature distribution. The outer surface tem- 
perature (boundary temperature) is known. Here the use of cylindrical 
coordinates facilitates analysis . In cylindrical coordinates the heat flow 



equation is: 

_L _I ( 1L\ + ±- * T .+ * 2 T = k A T (1) 
t k Mr J r 2 ^()> 2 ^z 2 At 

For the steady state , where = , then 

it 
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A 

^r 



}r 



i 2 T 



.2 ^2 



l 2 T 
172" =0 



(2) 



Since this deals with a nonhomogeneous cylinder, any approximation to 
the partial differential equations must show a dependence on the con- 
ductivity of the material. A suitable approximation may be derived by 
reference to the accompanying diagram. 
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First assume the boundary between materials to be midway between two 
points on the grid and introduce a boundary temperature T, . Then the 
quantity of heat which would pass from one square and (assuming incorrect 
initial assignment of temperatures) the quantity received by the second 
must be equal; or 

k i < T i " T b ) = k (T b - T > 



or 



T b = h T i + kT 
k + k,- 
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Then approximate 



or 



*r 1 



Vr\ 



= T b-T 
h/2 



TF i U+k 



Ti - T 

, h/2 



An equivalent approximation redefining T^ may be made from the opposite 
side to give: 



ir 2 



Jh + 1 



•i+1 + 



jH+i +k J \ h/2 / 

The final approximation to the first derivative may be taken as: 






iT 
*r 



i + (t) 2 MW- 



-T i+1+ T^ 



Vvrt 



The second partial may be taken as: 
* 2 T 



* r 2 



\^i\ \^<i \fer + k) (Stt) Vk)(S7r) 



For uniform boundary temperature: 
* 2 T 



and 



}0 2 
* 2 T 







iz- 



*l+3 +k 7 V h 2 /2 ) V k i+2 +k y Vh 2 /2 



T i+3 T \+ / k i+2 



T i + 2" T 



Making use of these approximations in equation (2) gives the steady-state 
heat flow difference equation: 



g +k )[ T 3- T ] + (g +k ) M 



= 



referring to the general cell illustrated. 
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Exercises 



With this equation the standard technique of relaxation by residuals may 
be used to find the steady-state temperature distribution. 

This problem has been discussed in order to illustrate how the relaxation 
technique may be applied in cases of two added types of complexities: 

1. Cylindrical coordinate system 

2. Nonhomogeneous material 



1. Write a FORTRAN program to perform a matrix multiplication. The 
maximum size of the input matrices is 20 by 20. The actual size of 
the matrix is L, m, n, and is to be read from data cards. 



Note: 



'ij 



A ik • % 



k=l 



*2. Using the method outlined in the matrix inversion topic, find the coef- 
ficient inverse and solution for: 

x 1 + 2x 2 + 3x 3 = 3 

2x-l + 3x 2 + 4x 3 = 5 

3x^ + 4x 2 + 4x 3 = 5 

*3. The discussion of the problem on the iterative solution for homogeneous 
equations mentions the use of the iterative process for nonhomogeneous 
equations . 

Given the set of simultaneous equations 

AX = B 
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one such iterative procedure is to rewrite the equations in the form 
X = CX + F 

and use 

X x = CX + F 
n+1 n 

as the iteration equation. 
For the set of equations 



25x 1 + 2x 2 + x 3 = 69 



2x 1 + lOx + x = 63 
x l + x 2 + ^3 = 43 



the iterative form is 



x x = -24XJ - 2x 2 - x 3 +69 
x 2 = -2xj - 9x 2 - x 3 +63 
Xo = -Xj - x 2 - 3xg +43 
Iterate for the values of Xp x 2 and x„ which satisfy the equations. 
*4. Perform relaxation on the below grid which has boundary conditions: 

XT rz "V" -f- V" 

boundary y 

for the potential described by Laplace's equation 

} 2 V . * 2 V 



+ 



°TC^ Ow2 



= 



Make use of symmetry wherever possible. 
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SECTION VII: EMPIRICAL RELATIONSHIPS 



The previous chapters have discussed various ways of handling problems 
adequately described by theoretical equations. Some attention has been 
given to each of the general steps in computer solution of the problems , 
namely: 

1 . Selection of the appropriate relationships . 

2. Mathematical manipulations involving algebra or the calculus to 
arrive at a suitable evaluation equation. 

3. Choice of numerical technique when necessary. 

4. FORTRAN programming and actual computer evaluation. 

It must be recognized that this approach does not suffice in situations 
where step 1 is impossible or where the known relationships do not 
adequately describe the problem. 

The lack of theoretical knowledge may often be adjusted for by the use of 
empirical data taken by means of experiments and test runs. 

The discussion of the vehicle simulation in a previous section introduced 
the use of tabular data and an empirically derived equation together with 
theoretical equations in the description of the vehicle operation. 

The engineer is no stranger to the use of tables of data. Page upon page 
of collected data arises in almost every test of an engineered part or 
system. Conclusions are as often drawn from the consideration of tables 
of data as from theoretical data. It is important that this source of infor- 
mation be available , as is , for use with computer analysis . At the same 
time the possibility often exists for derivation of an empirical relationship 
to fit the data. This section will consider some of the standard proce- 
dures for: 

1. Use of tables of data as they exist. 

2. Derivation of empirical relationships from the data itself. 

The discussion will be broken into four areas, as shown in the chart below: 





FUNCTIONS OF A 
SINGLE VARIABLE 


FUNCTIONS OF 
MULTIPLE VARIABLES 


TABLE 
LOOKUP 


a. Table storage 

b. Table interpolation 


a. Table storage 

b. Interpolation logic 


DATA 
FITTING 


a. Fitting Criteria 

1. Selected points 

2. Least squares 


a. Linear regression 

b. Multiple linear regression 

c. Nonlinear Estimation 

d. Dimensional Analysis 
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Functions of a Single Variable 



Given the following table of data, how might the proper value of Y for a 
given value of X be obtained during a computer run? 



X 



X(I) 



Y(I) 











0.0 


X(l) 


.913 


Y(l) 


4.0 


X(2) 


.930 


Y(2) 


8.0 


X(3) 


.941 


Y(3) 


12.0 


X(4) 


.946 


Y(4) 


16.0 


X(5) 


.948 


Y(5) 


20.0 


X(6) 


.950 


Y(6) 


24.0 


X(7) 


.951 


Y(7) 


28.0 


X(8) 


.948 


Y(8) 


32.0 


X(9) 


.944 


Y(9) 


36.0 


X(10) 


.938 


Y(10) 


40.0 


X(ll) 


.928 


Y(ll) 


44.0 


X(12) 


.914 


Y(12) 



TABLE LOOKUP 

One procedure is to load the entire table into the storage of the computer, 
and then for a given value of X search the table for the corresponding value 
of Y. 

A FORTRAN program to accomplish the loading of the above table into 
storage and table lookup for several values of X is shown in Figure 14. 
Note that the search is accomplished with the use of an IF statement 
within a DO loop. For the case of the argument X being exactly equal to 
a table entry value of X, the corresponding value of Y is simply selected 
and printed . When the argument X falls between two table entries linear 
interpolation is performed . The graphic derivation of the linear inter- 
polation equation is: 



fri+1 - yi) ( x " x i> 
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This is the equation of the straight line joining points (x^, y^) and (x^ + ^> yi+l) 1 
Evaluation of the right-hand side for a particular value of x gives the 
corresponding value of y. This equation is used in programming statement 4. 
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Figure 14. 

If the approximation of the function by a straight line in the interval 
tyi+l ~ Yj) is not of sufficient accuracy, the function may then be approx- 
imated by a parabola or higher-degree polynomial by such methods as the 
Lagrange interpolation formula. 

DATA FITTING 

If the problem involving use of this table were to be run many times on a 
computer, consideration may be given to finding an equation which will 
pass either through or within tolerance of all the points in the table. In 
most engineering problems it is sufficient to find an equation which passes 
within a specified tolerance of all the points in a table of data . 

Model Selection 

First a selection of the equation form to be fitted to the data must be made. 
A few of the possible selections are listed below. 



TYPE 



EQUATION 



1 . Polynomial 



2 . Logarithmic 

3 . Exponential 

4 . Power 

5. Fourier Series 



From y = slq + a.«x 

To y = a + aiXi + a. 9 x 9 + aoXo + a^x^ + a^x 



Generally restricted to this range, 
y = a Q + a 1 log x 

y 



l 5 A 5 



a„a. 



x 

"(Tl 

a l 
y = a Q x 



a„ + 




m 

Z 

n= 1 



(a cos nx + b sin nx) 

v n n ' 
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The actual decision may be made on the basis of theory, some preliminary 
plotting, past experience, or by trial and error. 

Model Fitting Criteria 

Next a criterion for actual fitting of the equation form (finding values for 
the a^) must be selected. Some of the possible methods are listed below: 

Selected Points . As many sets of observation data as there are aj to be 
determined are substituted into the selected equation, and the resulting 
system of equations is solved for the a^. While this method is very crude, 
it may be of value in situations in which available data is limited (see 
conical filter shell stress problem later in this section) . 

Harmonic Analysis . This widely used computer application is useful in 
fitting a Fourier series to a set of periodic data. 

Least Squares . This is the most widely used procedure for calculating 
the parameters a^ for the selected model. 

All of the models mentioned above (Fourier series excepted) may be fitted 
to a set of data by least squares . The principle will be explained with 
respect to the simple linear model 



y = a 



+ a-t x 



Given a table of n sets of data, 
where y is assumed to be the above 
linear function of x, determine aQ 
and a]_. 




The least squares criterion is to determine slq and a-^ such that the sum 
of the squares of the vertical distance between the data points and the 
straight line is a minimum. Referring to the graph, this may be stated 
mathematically as: 



L_, i = mi 



minimum 



i = 1 
This is true only it' 



» E £ i 2 

i = 1 



d a 



= 
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and n 



L* 



i = 1 
^ a 1 



= 



The sum may be expressed in terms of the equation to be fitted and the 
original data points . 

€i = Yi - (a + a x 34) 
Then: 



} L^ € i 2 ^ Z^ [y. - (a + a x Xj)] 

i = 1 i=1 =0 



3 a Q ^ a Q 



n n 



5 L € i 2 * L [>i-<a + a lXi )] 

i = 1 i= 1 

"} a^ > a^ 

Performing the differentiation and simplifying gives: 



11 " 

iL, *i - Lj a i x i - na = 

i = 1 i = 1 

n n n 

L^ x. y. - a ( /L 34 ) - a L x. =0 
i = 1 x i = l 7 i= 1 

These two linear, nonhomogeneous equations may be solved for a„ and a,. 
The parameters a.Q and a, , therefore, are computed in terms of sums and 
sums of cross products of the raw data. 



Functions of Multiple Variables 
TABLE LOOKUP 



Given a set of engineering data of the following form , from which x is to 
be computed for various sets of values of y and z: 
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x i yi 

*2 Y2 z l 

*3 Y3 

X4 Y4 

x 5 ys 

X6 Y6 z 2 

x 7 y 7 

x 8 Y8 

*9 y9 

x io yio 

xii yn z 3 

A table lookup with interpolation procedure , more complex logically to be 
sure, may again be used. This procedure may be stated as: 

(a) For zi interpolate for x as a function of y alone. 

(b) Store the resultant value of x along with z^ in another table. 

(c) Repeat steps (a) and (b) for all values of z giving a complete table of 
x as a function of z alone . 

(d) Interpolate in this resultant table for final value of x to be used. 

DATA FITTING 

Quite often, particularly in such fields as chemical engineering, the table 
lookup possibility is impractical for one of two reasons: 

1. The problem demands a fitting equation for a meaningful solution. 

2. The data cannot be obtained in a form similar to that shown in the 
previous table. 

When this situation is true, curve fitting may again be used. 

For example, an equation in the description of the vehicle simulation 
model was presented as an empirical relationship. In this case R is 
assumed to be a function of three variables M, A and V. 

2 
R = aM + bAV 

The least squares criteria may again be used in conjunction with exper- 
imental data to determine values for a and b which best fit the equation to 



the data. 

To apply least squares relationships in this instance would require calling 
AV 2 a new variable x 2 = AV , giving a linear relationship of the form: 

Y = ax- + bx„ 
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In general, data fitting of linear equations is called linear regression. 
If, as in this case, there are more than one independent variables, the 
procedure is called multiple linear regression. 

This chain of first assuming the form of a relationship and then using 
mathematical criteria to fit the relationship to experimental data can be 
thought of as a search for a useful "prediction equation. " 

A typical application of statistical techniques to arrive at experimental 
relationships may be represented by the problem of fuel oil smoking. 
Smoking is an undesirable trait of fuel oils . If the effect of smoking can 
be related quantitatively to properties of fuel oil, then smoking may be 
controlled. 

The method of attack is first to isolate the variables which affect smoking. 
In this case, they are aromatic content, olefin content, sulfur content, 
and boiling point. 

An arbitrary scale for measuring smoking experimentally must be agreed 
upon. Then, a large number of fuel oil samples are analyzed for the four 
characteristics above and burned to determine smoking. This gives rise 
to the following table of data: 

Boiling 
Sample Smoking %Aromatics %01efins %Sulfur Point 



No. 



X 1 



X, 



X, 



X, 



1 


4.3 


1.37 


.07 


.00 


42.9 


2 


3.7 


1.29 


.06 


.01 


41 


8 


3 


4.9 


1.41 


.13 


.02 


37 





4 


etc. 


• 


• 


• 







This, in turn, is subject to the following statistical analysis: 

1. Regression Equation: If smoking is assumed to be a linear function 
of the variables, the analysis will determine the best fitting equation: 

Y = bjX! + b 2 x 2 + b 3 x 3 + b 4 x 4 + b 5 
where bp b 2 , . . . , b 5 are the fitting constants. 

2. Partial Correlation Coefficients: The square of one of these coeffi- 
cients expresses the percentage of change in smoking which is due to 
a corresponding change in the respective variable. In this sense it is 
a measure of the dependence of smoking on the variable . 

3 . Multiple Correlation Coefficient: This gives a measure by which it 
may be determined whether all significant variables have been selected. 

4. Standard Error of Estimate: This gives the accuracy of prediction by 
the regression equation. 
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If the analysis indicates that the proper variables have been chosen and 
that the accuracy of the regression equation is sufficient, it should be 
possible to determine means to control smoking of the fuel oil. If not, 
other possible variables must be selected and new data obtained from 
which to determine other contributing factors . 

The preceding discussion has stated that the models chosen for multiple 
regression must be linear models. More explicitly, this means that the 
partial derivative of the model function with respect to one of the param- 
eters must be independent of that parameter. This mathematical 
restriction is severe for some desired applications. 

The Mathematics and Applications Department of IBM has developed an 
approach to fitting nonlinear models. The mathematical technique is 
interesting in that it consists of iterating through the two steps of (1) 
making a linear approximation to the nonlinear model and (2) fitting the 
linear model by multiple regression. The result is the possibility of 
fitting nonlinear models , an example of which might be 



= *1 



f a l x l a 2 X2 > 



a l + a 2 



This nonlinear estimation has been used further in fitting sets of equa- 
tions to data. 



Dimensional Analysis 



This technique for model selection of relationships between physical 
variables has wide application to engineering problems . With electronic 
computation the application of dimensional analysis can be greatly facil- 
itated. 

As usually written, each physical variable of an engineering equation is 
expressed in some combination of fundamental units. These units — time, 
length, mass, temperature, etc. — are the dimensions. 

Dimensional analysis is a systematic study of the dimensions of an equa- 
tion containing several variables in order to: 

1. Discover any variables whose dimensions are incompatible. 

2. Ascertain the nature (but not the true value) of the relationship between 
the variables. 

3. Assemble the variables in dimensionless groups so that the relation- 



ship between these groups can be expressed algebraically. 

As an engineering application of dimensional analysis which will make the 
procedure clear, consider a problem of conical filter shell stress. A 
filter shell must be designed to withstand large stresses applied through 
a tightened nut at the top. This force is needed to overcome the force due 
to oil under pressure inside the shell and maintain a firm seal at the base 
gasket. 
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ZZA 



IHJ 



einforcing 




Bolt 



IHJ 



Shell 



Base 
Gasket 



An equation is desired which will predict the applied stress, W, for which 
the shell will fail. For this structure a complete theoretical analysis is 
too involved. Therefore, an empirical relationship will be determined. 
It is assumed that variables involved are: 

S = Ultimate shear stress of the material 

t = Thickness of the shell 

D = Mean diameter of shell 

■&■ = Cone angle with the vertical 

d = Diameter of reinforcing ring 

W = Applied stress 

Dimensional analysis now consists of finding dimensionless groups of 
variables so that a dimensionally correct equation may be written. An 
approach here is to let 

tt x = S a t b D c d 



it 2 = S e t f D g -e- 



7T 3 = S n t 1 D) W 



(1) 



The dimensions of each of the variables to be inserted are: 



S = MT" 2 L" 1 
t = L 
D = L 
d =L 

■e- =1 

W = MT -2 L" 1 



M for mass 
L for length 
T for time 
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For tt i to be dimensionless, the sum of the exponents for each funda- 
mental dimension must equal zero . 

For M, a = 

For T, a = 

For L, -la + lb + lc + 1 = or b = -(1+c) 

For this particular case the resulting three equations have an infinite 
number of solutions. For simplicity consider: 

c = 
Then 

b = -1 

and -rr 1 = S t _1 D° d 

Or 7T , = (A.) 
1 t 

Similarly, w 2 = (Jl) -0- 

(Since -9- is dimensionless, this ir factor may be rewritten as two dimen- 
sionless groups.) And 

3 w 

At this point the variables are arranged in dimensionless groups and a 
relationship between the groups may be assumed. 

A reasonable form to assume is: 



3 - v ■• 2 ) a i ■■ ! 
or 



7T.=K(7T )- (7T V 



St 2 _^/D\a /_d_\b .©? (2) 



w -*(fr (f-r* 



where -0- has been used as a separate tt factor because it is dimension- 
less. 



Rewriting (2) gives: 

g = KW_ /D_\a /d\b ^c 



t 2 U / \t / (3) 

where K, a, b and c are constants to be determined. 

One method for calculating suitable values for the constants K, a, b and 
c is that of selected points covered earlier in this section. 

(a) Take the log of each side of (1) 

KW n a 

log S = log (-£-) + a log (■—) + b log (-S-) + c log -6- (4) 
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Probabilities 



(b) Make four shells of different dimensions with known ultimate stress, S. 

(c) Crush shells to obtain W. 

(d) Substitute all the values into (4) to obtain four equations in four un- 
knowns (K, a, b, c) and solve these equations simultaneously. 

(e) Apply constants to (3) and verify by testing other shells . 

For improved accuracy at the expense of shells , the following procedure 
may be used (this is equivalent to the least squares done graphically): 

(a) Holding S, t, D and d constant, crush many shells of varying -Q- . 
For this case (4) becomes 

- log W = c log -0- + B * 
where B is a constant. 

(b) Plot W versus -9- on log-log paper. The slope of the straight line is 
then c. 

(c) Repeat a and b for other variables. 

Note the possibilities for computer usage on larger problems of this kind. 
Simultaneous equations (to be discussed later) are handled twice. Further- 
more, the procedure for improving the accuracy is a simple curve fitting 
problem. In fact, the whole procedure for dimensional analysis, since 
it is well defined, can be programmed for a computer. 



The probability of a variable having a certain value is a question which 
often arises in engineering. If historical data on the variable is available, 
it is often possible to do the necessary predicting. This problem illus- 
trates a method for handling normally distributed variables. 

Given a table of observed values for a particular variable whose frequency 
distribution is assumed to be normal: 

(a) Find the probability of a reading being less than any given value of 
x, and 

(b) Find the probability of a particular reading lying between any two 
values of x. 
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The general distribution function is: 



1 /x-m \ 2 



f(x) 



2\ a- 



277- 



where m = mean 

ar = standard deviation 



The standard normal curve is given by: 

e 
f(t) = = 

V2 7T 



-t 2 /2 



Therefore a transformation of variables will allow us to make use of 
standard normal distribution function. This transformation is 



t = 



x-m 



The probability of a reading falling below a given value of x is given by: 



.x 



J f(x)dx=Jf(t)dt 
-00 



-03 



x-m 



where t = 

The probability of a reading lying between x^ and x 2 is: 



J f (x)dx - f f (t)dt 



x n -m 

where t.. = _ ± 

1 cr 



t = x 2~ m 
2 ^ 
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Exercises 



The problem, then, is one of evaluating the integral 



b ,2 



9-rr J e 



V2= > e " dt 



a 



A table of values of this integral for various values of t (normal areas and 
ordinates) is available in any statistics text. In computing, problems of 
this type might be tackled in two ways: 

1. Store the table in memory. 

2. Use numerical approximation for evaluating the integral. 

An approximation by Hastings* might be used with the latter approach. It 
is 

x 

2 



jzjf J e dt for o < X < 00 



~ 2 f -t 
1 



where (J) (x) = 1 - 4 

[l+a-jX+agX +a 3 x 3 +a 4 x J 



and 



a = .278393 a = .000972 

a n = .230389 a. - .078108 

2 4 



1. Write a FORTRAN-language program to perform table lookup with 
linear interpolation for a tabulated function of two variables as 
described under "Functions of Multiple Variables." 

2. Evaluate the coefficients for the linear least squares fit 

y = a Q + a^ 
for the following table of data 



* Approximations for Digital Computers by Cecil Hastings, Jr. , 
Princeton University Press, 1955. 
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X 



14.0 64. 

13.9 54. 

13.5 55. 
13.4 56. 

12.6 46. 
12.6 51. 
11.8 42. 

11.4 47. 
11.3 48. 

10.5 35. 
10.3 44. 
10.3 29. 

9.3 31. 

7.6 14. 



*3. The distribution in measured diameter of the dirt particles normally 
deposited in a carburetor sediment bowl is a normal one. 

For a mean of .0025 centimeters and a standard deviation of .0002, find 
the probability of a dirt particle being greater than . 0027 centimeters in 
diameter. 

Use any statistics text for necessary tables. 
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SECTION VIII: OPERATIONS RESEARCH 



The reason for including this as a topic of engineering analysis for com- 
puters is that often the technical personnel of a company will be turned to 
for answers concerning the mathematical aspects of Operations Research. 
Books have been written on each of the subjects to be mentioned. The 
intention here is merely to present examples in such a way as to make 
possible the recognition of problem types as they occur . 

Operations Research can be thought of as the application of scientific 
methods, techniques or tools in the management decision area. It attempts 
to present a clear picture of the system under study and at the same time 
predict the consequences of alternate sets of actions — in other words, 
show the best decision. Methods for attaining either or both of these 
goals require knowledge of the system . Managed systems are , in general , 
complex, which means that vast amounts of data are required to describe 
these systems . Manipulation of this data to attain the stated goals gen- 
erally requires a computer. 

OR seems to be at the "clear picture" stage at the present time. The more 
complex problem of "prediction" has been hindered by the fact that 

1. The mathematical and logical methods employed by OR are unfamiliar 
to those most intimately concerned with the problems (the executives). 

2. The historical data describing the system does not exist in proper or 
even intelligible form. 

However, consideration of the prediction level of OR while working on the 
clear-picture level has the following benefits: 

1. Creating a clear picture of a system results in the organization of 
data into intelligible form. This data can later be used at the pre- 
diction level of analysis . 

2. Consideration, even without use, of more advanced techniques often 
gives a better understanding of, or "feeling for," a system. 

Almost any procedure for presenting data to an executive can be thought 
of as part of the clear-picture operation. However, the computer allows 
concise data handling in many new areas at reasonable costs. Some of 
these are: 

1. Graphic presentation of parts failures in the field, showing effect of 
engineering changes on failures, failures per unit of operation, fail- 
ures as a function of production group, etc. , on a day-to-day basis. 

2. Presentation of historical data concerning use of spare parts in the 
past so as to gain some insight into what an "all-time production" of 
a spare part should be. 
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3. Creation of a weekly report to a supervisor, showing jobs assigned to 
his shop, estimated time for each job, group within shop to which it was 
assigned, actual time for completed jobs, and year-to-date figures for 
number of jobs assigned to each group with totals for estimated and actual 
job times. This would allow analysis of each group's performance. 

What are some of the prediction techniques presently being employed? 

1. Linear Programming: This technique may be used on any system 
which may be represented by a set of linear equations. 

There may be more variables than equations and therefore an infinite 
number of solutions to the set of equations. The object of linear pro- 
gramming is to find a particular solution which will optimize (make 
maximum or minimum) another linear function (usually a profit or 
cost function). This technique has been very successful in the blend- 
ing of petroleum products and the mixing of grain feeds . It is not un- 
usual for the mathematical model of a gasoline refinery to be as large 
as 100 equations in 500 variables. The actual mathematical techniques 
employed for linear programming will not be discussed here. However, 
it is worth noting that the techniques have been refined to the point 
where their application can be rather mechanical. That is to say, a 
knowledge of the mathematics is not required for successful use of 
linear programming. 

An important special case of linear programming is the transporta- 
tion problem. It may be used to optimize any system for which a com- 
modity is available at several sources and is required at several 
destinations , and for which the transportation rate for each possible 
source-to-destination combination is known. The most economical 
shipping assignment will be determined. 

2. Forecasting: Historical data is analyzed statistically to enable simul- 
taneous forecasting of many interrelated factors of management inter- 
est in terms of other known variables that affect them. Although, in 
general, the statistical techniques are not new, the computer allows 
consideration of many more variables and, therefore, the handling of 
complex management systems . 

3. Simulation: Here a mathematical and logical model of the system in 
question is created. This model is operated according to the normal 
rules of the system operation for a specified period of time. Con- 
cise data on the results is made available. This gives a picture of 
how the system will operate — quickly, and in clear, concise form. 
More important, the system may be altered or the rules changed and 

the model rim again; thus, the effect of alternative decisions can be 



tested on paper. 

The term "simulation" was introduced earlier in the discussion of 
vehicle simulation. In the present context, however, the system 
simulated can be a company organization and/or its manufacturing 
tools — • a much larger concept. 
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Some discussion of the techniques of simulation is in order here, since, 
while the problems tackled via simulation may be extremely complex, the 
basic mathematical elements are quite simple. 

There are two basic elements of any simulation, the first being the repro- 
duction of actual (or perhaps expected) events . This can be done either 
by using actual historical data (for example, past sales record by month) 
or by using statistical parameters which describe the desired behavior. 



To reproduce arrival of orders by size when the past 



Frequency 
Distribution Curve 




frequency distribution of order sizes is known, an accumulative distri- 
bution curve is constructed as the first step. 



100 



% of orders 
greater in 
size 



Accumulative 
Distribution Curve 




Using random numbers (generated or taken from a table) to represent the 
M % of orders greater in size," the corresponding order size can be read 
from the curve or obtained from a table or equation representing the curve, 

This use of random numbers and the accumulative distribution curve will 
reproduce the behavior pattern from which the curve was constructed. 
Further, each order size will occur at random within the overall pattern. 
This is exactly what is desired in simulation. 

The second element in simulation of a system is time . Two distinct 
approaches have been used here: 
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1 . The system is reviewed at definite intervals of time and all activities 
are updated at each review. A disadvantage of this approach is that 
for intervals of review, no activity may have occurred. However, 
many systems with natural periods can be simulated nicely using this 
procedure . Simulation of inventory and truck fleet operations fits this 
pattern. 

2 . The system is reviewed only when something important happens . These 
important points in time may be recorded in a time status record 
(TSR) . Initially a distribution-of -events curve is consulted for the 
entry to the TSR. Subsequently each time an entry to the TSR is han- 
dled by updating the activity, the first step is to sample the events, 
distribution for the next TSR entry. 



Job Shop Simulation 



A job shop is distinguished from an assembly-line shop by its flexibility 
in processing articles of manufacture. The amount of time required for 
a given operation varies from article to article and the route of operations 
varies from article to article. 

Important advantages may often be realized through simulation of a job 
shop on a computer, considering such things as machines available, labor 
available, operating costs, fluctuation in order s , priority of orders , 
operating rules, etc. Records of order completions, in-process inven- 
tory, and waiting lines in front of machines can be maintained. The effect 
of proposed changes in the physical plant or operating rules may be pre- 
tested using such a simulator. 

A short description of the logic which goes into the simulation of a single 
operation job shop will illustrate the basics of this type of simulator. It 
illustrates the use of probability distributions and the time status record 
approach of handling the time element. 

The logic diagram for a single operation job shop might be as follows: 



Finished 



Pieces Arrive 






The information required to simulate this system would be: 

1. The probability distribution of time between the arrival of pieces. 



2. The probability distribution of operation time per piece. 
The results desired would be: 

1. The number of pieces which can be completed per unit time. 

2 . The amount of idle time experienced per unit time . 
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3. The average waiting time per piece. 

On the basis of the two probability distributions, a time status record is 
generated which guides the simulation program. The procedure for ma- 
nipulating the model is shown in Figure 15. 







Piece is 
completed 
in A 



Move piece 
from A to F 




Move piece 
from W to A 



Place piece 
in A 



Sample for 
operation time 
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inTSR 



START 





Piece 
arrives for 
service 



Sample for 
next arrival 
and record 
in TSR 




Place piece 
in W 




Operation 
Completion 



Arrival 
of piece 



Figure 15. Job Shop Simulation 
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Two curves which might represent the required distributions are shown 
in Figure 16 along with a short, random sequence of numbers. Figure 17 
shows the time status record, and the associated recordkeeping which 
might be performed by the computer in simulating the operation. 



The conditions for beginning this simulation are: 

(a) Line 1 represents the status of the one operation shop at the beginning 
of the day. 

(b) Lines 2 and 3 are entries to the TSR made during the previous day. 

(c) Lines 4 and beyond represent entries generated during this day's 
simulation. 

(d) The simulation of this day's events begins at 9:00 a.m. (the first event) 
and can be followed in the logic chart beginning at the point marked 
START. 

This simple statistical technique for accomplishing simulation can be used 
for a wide variety of operations research problems. 
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Figure 16 
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TIME STATUS RECORD 
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Inventory Control 



Figure 17 



Inventory control is important to industry for sizable dollar reasons . It 
is an important area of effort computer-wise because (1) much of the 
recordkeeping required for inventory control is presently being done by 
data processing systems, and (2) many problems of inventory control can 
be tackled by those mathematical treatments which are greatly facilitated 
by use of computers . 

Inventory control, from a data processing viewpoint, is best broken into 
two distinct procedures: 

1. Forecasting the future sale of items from inventory. 

2. Use of the resulting forecast figures in a control procedure which will 
maintain that level of inventory which best meets the requirements laid 
down by management operating decisions . 

To take a hypothetical company in order to illustrate a possible study of 
these two areas, assume that the Best Parts Company has a central ware- 
house which supplies retail dealers . The situation contains the following 
factors: 

1. The inventory status must be reviewed each week. Total sales for each 
item are recorded and needed replenishments of stock are ordered at 
each review. 

2. In practice, all replenishment orders to the plant can be filled within 
a lead time of 1 1/2 weeks . 

3. The activity of different parts varies greatly. 

4. There are no largely seasonal parts. 
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5. Two years of sales records are available. This record gives the 
total sales of each part for each one-week period. 

Best Parts Company management believes that it can improve the control 
of the inventory in one or more of three ways: (1) reducing inventory, 
(2) reducing out-of-stock conditions, and (3) giving simpler processing 
of control information. 

A study of the control procedure is initiated with the understanding that 
the computer facility, with whatever new techniques available, may be 
used. The study takes the following form: 

Forecasting 

Four approaches to arriving at the forecast sales to retailers for each 
part are considered: 

1. The use of linear regression in finding a prediction relationship of the 
form 

Forecast Sale = a + a^x^ + a£X£ + a-i^n 

where the x's may be factors felt to correlate with or influence the 
sales, such as gross national product, number of people in $5,000 - 
$10,000 income bracket, etc. This approach requires extensive 
analytic effort, much diverse data, and large amounts of computing 
effort. 

2. Classical time series analysis gives a forecast based only on analysis 
of previous sales data. The linear trend line of the form 

Forecast Sale = a + bT 

where T is time in sales periods, requires only past sales data, but 
considerable computing for fitting of a and b. 

3 . A common practice in industry is to use a moving average for sales 
forecasting. The formula for a four-week moving average is 

F „=S+S ., + S n + S „ 
r+1 r r-1 r-2 r-3 



where F , n is the forecast sales for week r+1 and the S's are the 

r+l 
actual sales for the previous four weeks . This method is simple , 

data-wise and computation- wise. A major disadvantage is that it does 

not generate statistics which are found important in the control pro- 

cedure — in particular, the statistic called the standard deviation of 

actual sales from the forecast. Approaches 1 and 2 may generate 

this statistic. 

4. A weighting technique credited to Robert K. Brown and called exponen- 
tial smoothing has gained attention and application in the past few years. 
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Forecasting is based, as in the moving average, on past sales data 
and requires even less computing. A major advantage is that a sim- 
ply computed estimate of the required standard deviation is available 
and corrected to the most recent sales review figures. 

The exponential smoothing procedure is chosen for study on the basis of 
limited available data and limited computing time. The forecast equations 
for exponential smoothing may be given as 

P i = P i-1 +oC < S i-l " P i*-1> (Demand) 

Tj = ^ (Pj - Pj*.!) + (1 - oC ) T.^ (Trend) 
P i* = P i + ^ZF^ < T i) (Forecast) 

D i = °< (S i-l " P i*-1 ) + (1 " * > D i-1 (Deviation) 
where 



Then 



and 



with 



S^ = sales figure for period just ending. 



Pi = demand expected for next period 

Ti = trend measure of whether last several periods have seen a 
rise or fall 



Pi* = trend adjusted, final forecast, for next period sales 



1.25 Di = estimate of standard deviation of actual sales from fore- 
cast. 

The constant o< represents the weighting Choice. Without the trend ad- 
justment, the demand (Pi) computed will closely approximate the moving 
average based on N periods if oC is chosen as 

2 

With the use of the trend adjustment, experience has shown that a choice 
of «< = . 1 is most satisfactory. 

Control Procedure 

The key to inventory control is in the way one parameter is computed. 
This parameter is the requirements figure. The requirements figure 
must be arrived at in such a way that it relects (1) the knowledge of actual 
sales, so that stock is available to meet the demand, and (2) management's 
ideas of what is adequate by computing the proper stock level, within a 
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defined approximation, so that out-of-stock conditions are kept below a 
set level. 

This implies that management must set the level for out-of-stock condi- 
tions — and this is true. The familiar quality decisions such as "we must 
have few out-of-stock conditions" or "we must reduce inventory" must be 
replaced with decisions such as "an average of 5% out-of-stock is to be 
our goal." 

Literature is available on methods of arriving at a proper out-of-stock 
percentage figure which tends to minimize inventory operating costs for 
particular inventory situations. 

In this example, the management percentage level figure is assumed 
available . 

One approach to defining the necessary requirements equation is to con- 
struct an inventory model similar to that shown below. This model assumes 
the following "ideal" conditions: 

1. It represents a possible inventory level as a function of time for one 
item. 

2. The sales are linear — the same sales figure for each period. 

3. The sales in each period are exactly the amount forecast for that 
period. 

4. The lead time is 1 1/2 review periods (lead time being the time be- 
tween generation of an order and actual receipt of replenishment 
stock) . 

The graph begins at a point in time where the actual inventory is 
exactly at the requirements level. 



Requirements Level 




Time 
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The upper curve represents the inventory actually in the warehouse 
at any given time plus the inventory on order. 

The graph shows that for a requirements figure computed strictly in terms 
of the forecast and the lead time, zero stock conditions are to be expected. 
Since the Best Parts Company knows that it is not operating under "ideal" 
conditions (sales fluctuate from period to period), a safety stock is needed. 
The final requirements equation then is: 

R i = p i* + Lead * p i* + Safety 
The Safety Stock Calculation 

Reference must be made here to the discussion of probabilities in the 
chapter on empirical relationships. If a frequency distribution of errors 
in the forecast over a period of time is plotted, it might look like Figure 
18 (assuming a normal distribution) . 

The accumulative area under the curve, as we move from left to right 
along the error axis, represents the probability of the error in forecast 
being equal to or less than the corresponding value of the error. 




-l o 

Error 
(Forecast — Actual Sale) 

Figure 18 

If the shaded area in Figure 18 represents 5% of the total area, then only 
5% of the time is the error in forecast greater than four units. 

Therefore, if four units are always added to the inventory stock require- 
ments, the out-of-stock conditions should occur only for errors larger 
than four units — or only 5% of the time . 

For computational purposes, the same proper error value for a stated 
probability may be obtained by using the transformation equation: 



t = 



x-m 



where for this use 
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Then 



t = value corresponding to required probability level and is taken 
from the tabulated values available for the normal distribution 

x = desired error value 

cr = standard deviation of error 

m = mean of the errors (in this case m = o) 

x=crt 



and since 1.25 D^ is an estimate of cr, the safety stock required may 
be computed as 

x = 1.2.5 Dj *t 

The value of t corresponding to a probability level allowing 5% out-of-stocks 
is 1.65. The final requirements equation to be used by Best Parts Com- 
pany is 

R i = p i* + Lead * P i* + 1 - 65 " 1 - 25 ' D i 

Note that the statistic D, like all the other quantities here, is computed 
for each item in the inventory. The requirements figure for an item level 
is thus dependent on that item's own past sales behavior. Items with steady 
sales will have small safety stocks. Items with widely fluctuating sales will 
have large safety stocks. This may often make possible the seeming para- 
dox of reducing overall inventory and at the same time reducing out-of-stock 
situations . 

This completes the two basic selections needed in establishing an inven- 
tory control procedure — a forecasting technique and a requirements com- 
puting technique . All other quantities relating to inventory control can be 
computed each period in terms of the forecast and the requirements. For 
example, the replenishments equation for the inventory might be given as: 

ORDER = REQUIREMENTS - ON HAND - ON ORDER 

This procedure requires extensive testing before it may be successfully 
initiated. One way of accomplishing this testing is to simulate on a com- 
puter what would have happened if the procedure had been used over some 
time interval in the past. A simulation of inventory control has two basic 
elements: actual events and time. Actual sales data (quantity of sales 
of each item or each period in the time interval) is all that is needed to 



satisfy the reproduction of actual events elements. The normal use of a 
review period in inventory control makes it a natural for the use of re- 
view at intervals in handling the time element. 

A logic flow chart for a computer program to simulate the inventory model 
for a single inventory item is shown in Figure 19. 
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An important result of such a simulation is the data needed to prepare 
charts similar to that shown in Figure 20 . 
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Figure 20 



A graphic display of this kind, showing what will happen if a procedure is 
adopted, represents an excellent pretesting of the practicality of the pro- 
cedure . 



Linear Programming 



Mathematically, linear programming is finding a solution of a set of 
equations similar to: 

a ll x l + a 12 x 2 + a 13 x 3 + a 14 x 4 + a 15X5 = b l 

a 21 x l + a 22 x 2 + a 23 x 3 + a 24 x 4 + a 25 x 5 = b 2 

a 31 x l + a 32 x 2 + a 33 X 3 + a 34 x 4 + a 35 X 5 = b 3 

in which we have more variables than equations, and therefore an infinite 
number of solutions . 

The particular solution chosen is the one which maximizes or minimizes 
an additional equation (generally a cost equation) which might look like: 

cost = CiX-i + CgXo + c q x 3 + C 4 X 4 + C 5 X 5 



Any system (either physical, such as a petroleum refinery, or operational, 
such as work assignments, or a combination of both) which may be de- 
scribed in terms of linear relationships may be studied by linear program- 
ming. 

Many relationships which are not linear may be approximated by combi- 
nations of linear segments as illustrated graphically in the accompanying 
figure. 
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f(x) 




Much of the present-day power of linear programming is due to its having 
been taken out of the hands of the mathematicians and reduced to a fairly 
common engineering procedure. Rather than working with sets of equa- 
tions, the engineer works with an array or matrix which is the format re- 
quired for entering data to the computer. In place of mathematical opera- 
tions in creating this matrix, the engineer employs straightforward rules. 
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1.11 


1.25 


-.15 


-.2 










.5 


Reformer Capacity 










1.11 


1.25 


.05 












.6 


Butane 








.1 


-.03 


-.04 




-.01 










.05 


Straight Run 
Gasoline 


1.0 












-.1 


-.15 










.15 



Figure 21. 



A look at the array in Figure 21, which represents a complex problem to 
be optimized, will illustrate the use of such rules. 

This array represents the statement of the problem to be optimized. The 
objective will be to determine a recipe for blending available petroleum 
components to obtain a motor fuel meeting several specifications — and 
return maximum profit to the refinery. 

Each horizontal row in the array represents a restriction (equation) — and 
each vertical column a variable. The actual numbers in the array are 
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coefficients which would appear in the corresponding set of linear equa- 
tions which are required to describe the problem mathematically. The 
variables are the components to be blended. The restrictions may be of 
two types: 

1. Quantity restrictions — for example, where the amount of motor fuel 
to be produced is restricted. 

2. Qualitative restrictions — for example, where the Reid vapor pressure 
of the resulting fuel is restricted. 

Line 4 in the body of the array is read: 1.0 times amount of straight run 
gasoline to be used in the blend, plus 1.0 times amount of alkylate to be 
used in the blend, plus — — ~ — - — •, minus 1.0 times excess motor fuel, 
must be equal to the 1.0 barrels of motor fuel which is required. The 
coefficients can be entered to the array without ever putting down the 
equation stated above for this quantitative restriction. 

Row 1 represents the representation of the octane number restriction. 
The refinery sets a minimum level for the resultant octane number. In 
this case, let the minimum specification be S and let the actual octane 

ratings of the variables Xj to be blended be Ij for i = 1, 2 n 

where n is the number of variables. 

The assumption is that the resultant blend octane rating will be a weighted 
average by volume of the individual specification or: 

n 
S + P = I Ii Xi 
i=l 



n 

i=l 

where the P allows for the weighted average to be greater than the mini- 
mum specifications. A suitable form of the above equation is: 



Xi = 



n 
The term P X Xj is called the slack. 
i=l 

Note that the coefficients for each of the variables is to be S-Ij or the 



n 




n 


i 


(S - Ii) Xi + P 


i 


i=l 




i=l 



minimum specification minus the individual specification. 

This characteristic allows two rules for automatically entering coefficients 
into an array to handle most quality restrictions: 

1 . The magnitude of the coefficient is to be the absolute value of the 

difference between the required specification and the individual speci- 
fications . 
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2 . The sign of the coefficient is chosen positive if the individual specifi- 
cation exceeds the required specification, negative if the individual 
specification is below the required specification. 

Linear programming codes exist for most of the digital computers in 
commercial use. Once the array has been determined to adequately de- 
scribe the system under study, the data represented by the array is pre- 
sented to a computer along with the linear program. 

The output to be expected from the optimization computation is: 

1. Variable activities. This is a list of those variables chosen to satisfy 
the restrictions . With each variable is given the amount of the variable 
to be used in the resultant product. 

Many linear programs also indicate the range, if any, through which 
this amount may be varied, leaving the solution at a maximum or 
minimum as the case may be . 

2. Shadow prices. For those variables not selected as part of the solu- 
tion, an indication of the reason is given by the shadow prices. The 
shadow price for an unselected variable represents the cost of forcing 
one unit of the variable into the solution, considering the adjustments 
that would be necessary on other variables to again satisfy the stated 
restrictions . Again there is often indicated a range through which 
the cost or profit coefficient for the unselected variable may be moved 
without removing the optimization. 

Results of linear programming computation are used as a recipe for 
day-to-day selection of alternatives in such diverse areas as gasoline 
blending of available stocks, mixture of grain for feed, blending of 
meats for sausage, and slitting of paper rolls and fabric stock. 

The technique is also used on a study basis to aid long-range planning 
in the selection of refinery crude stocks , the selection of machine 
components in process industries, the scheduling of machine loads, 
and other areas. 

SUMMARY ON OPERATIONS RESEARCH 

As a brief summary of Operations Research, there are three techniques 
which have gained considerable acceptance and use in the study of man/ 
machine systems: linear programming, forecasting, and simulation. 

Linear programming computes the optimum choice of decision alterna- 
tives for a linear system within the available definition of that system. 

Forecasting by statistical means has important uses in long-range plan- 
ning and in day-to-day control of such things as inventory. 

Simulation may be used to study systems thus far undescribed mathemat- 
ically, by reproducing the logical operation of the system in testing alter- 
nate decisions . 
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Exercise 



1. A company would like to simulate the operation of a fleet of trucks 
which deliver from one factory to many parts of the country. The 
delivery is to be made for one product. The company has the follow- 
ing data which is important to the simulation: 



Day 



Number of Orders 
(Units of 1 Truckload) 



150 



100 
97 

115 
40 
80 



102 



Number of Days 
Required to Deliver 
and Return 



Number of Times 

This Period Was 

Required 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 



200 

500 

700 

2000 

4000 

5500 

3500 

3000 

500 

300 

40 





a. Sketch the distribution curve which might be used in simulation of 
the receipt of orders. 

b. Sketch the distribution curve which might be used to determine the 
length of time required to deliver each order (turnaround time) 
during the simulation. ' ^__ - 



The company would like to determine the best number of trucks to use 
to satisfy the delivery requirements while keeping the truck fleet cost 
as low as possible. 

The company estimates that it costs D dollars to maintain and operate 
one truck each year, and that the cost of back orders in terms of future 
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sales lost is: 

C = ANi + BN 2 

where A and B are constants, N^ is the number of back orders not 
more than seven days old, and N 2 is the number of back orders over 
seven days . 

Deliveries may be made every day of the year. 

c. Draw a logic flow chart of a simulation program to operate for one 
year, which the company might use to determine the optimum num- 
ber of trucks to maintain in the truck fleet. 
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EXERCISE 6 

Generalized flow chart to compute F(x) = aj x n + a2X n " 1 +a n +i 





READ 










N, A. 






I = 1,2- 


— ,N+1 






ir u 






READ 






X 






ir 






Set 






1= 1 






F(x) 


= 0. 






\ 


r i 


r 










Compute 
F(x) = F(x)-X + A. 




Step 
1 = 1+1 




i 


I 




\ 


f 






./Test ^v. 


1 No 1 




S. I = N+ 1 ?^^ 


K2J 




1 


r 


i 


ss ) 




PRINT 






X, F (x) 
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Section V 



EXERCISE 3 

The initial conditions may be used to complete the solution at t = . 



*< = o 



dx^ 
dt 



B -^n-(f 



(S) o — ^o-^o(|) 



Succeeding values for the variables and their derivatives may be obtained 
by two-level iteration. For each value of t, iteration is used to find a 
solution to oC = t + sin oC t. The following logic diagram makes this clear: 





0. 1, 2, 3, 4 





Compute 






y n = y n 








(-) " 


eu 


v dt 2 


u 


At 


£l- 


B -«< - 








#v 


A-xy 


--( 


-) 








dt' 
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Section VI 






EXERCISE 2 




Inverse is 




-4 4 -1 




4-5 2 




-1 2 -1 




EXERCISE 3 




x 1 = 2 




x 2 = 5 




x 3 = 9 




EXERCISE 4 



Solution is 



x x = 3 



x 9 =-3 



Xo = 2 



/28 


/28 


y28 


/31 


/37 


50 

41 


33 


34 


29 


30 


29 


28 


28 


28 


26 


28 


28 


28 


27 


25 
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Section VII 

EXERCISE 3 



t = x - m _ . 0027 - .0025 = i.q 
a .0002 

for t = 1.0, from table of probabilities 
CD 



-| V 27T 

I 



2.5 



2 
e dt = . 1587 



f -tV2 
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The Data Processing Bibliography (form J20-8014) provides abstracts of 
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